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ABSTRACT 


Near-vertical  incidence  deep  reflecting  events  from 
within  the  earth’s  crust  have  been  recorded  on  seismograms 
along  three  different  lines  in  southern  Alberta.  Predesigned 
systems  of  geophone  and  hole  patterns  have  been  synthesized 
into  an  effective  filter  for  long-period  surface  waves  while 
vertical-  and  horizontal-component  geophones  established 
their  vertical  incidence  character. 

The  data  were  recorded  on  FM  magnetic  tape  recorders 
with  a  wide  dynamic  range  and  the  results  have  been  digitized. 
By  application  of  digital  processing  techniques,  the  signal- 
to-noise  ratio  has  been  enhanced.  Power  spectra  calculations 
indicate  that  the  reflected  wavelets  have  energy  contained 
in  the  range  5  to  15  cps . 

A  synthetic  seismogram  program  has  been  developed. 
Three  important  results  have  been  obtained  from  its  applica¬ 
tion:  (l)  multiple  reflections  are  severely  attenuated 

within  1.5  seconds  after  a  primary  event;  (2)  any  transition 
layer  yielding  a  deep  reflection  must  be  less  than  1  km.  in 
extent;  (3)  the  Conrad  velocity  discontinuity  is  a  transition 
over  less  than  0.5  km. 


M  '  •  n .  't  ■  9  r  ^nol 

.  rtol  svi;l  9  n.G  olnX 

i  M1'  .  •  •• .  •  •  ; 

!  ;o  *  -our 

i  sri  :  <■  £ i j  .  ic'  o  ..«■  .  v  ;T 

J  '  0  0lr'  ■  /  J  ’  :j  1 '  .  •  * 

•  flltf  £  •  r  >’  ■'  ;  r’  " 


iv. 


A  strong  reflection  at  11.6  seconds  representing 
the  Conrad  discontinuity  has  been  continuously  correlated 
over  30  km.,  enabling  average  vertical  velocity  and  depth 
values  to  be  calculated.  Other  reflections  are  evident  on 
some  seismograms  suggesting  that  the  deeper  crust  may  be 
more  complex  than  can  be  resolved  by  refraction  methods. 

In  addition,  evidence  is  presented  which  suggests  a  high 
angle  fault  deep  within  the  crust  of  the  earth. 

It  is  hoped  that  this  thesis  will  demonstrate  the 
applicability  of  the  seismic  reflection  technique  for  deep 
crustal  research. 
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CHAPTER  1. 


INTRODUCTION 


1.1  The  Crust 

Knowledge  of  deep  structure  within  the  crust  of 
the  earth  would  greatly  enhance  our  understanding  of  its 
historical  evolution.  With  such  knowledge,  it  might  even  be 
possible  to  predict  the  occurrence  of  some  earthquakes  and 
so  minimize  the  destruction  they  cause.  The  seismic  reflec¬ 
tion  method  of  geophysical  research  is  particularly  suited 
for  yielding  structural  information,  as  evidenced  by  its 
widespread  use  in  petroleum  exploration.  However,  its  use 
in  exploring  deep  crustal  structure  has  been  limited  by  tech¬ 
nological  difficulties,  but  these  are  now  being  overcome. 

This  thesis  examines  various  methods  by  which  reflected  energy 
may  be  enhanced  over  different  types  of  disturbing  noise  and 
gives  details  about  crustal  properties  obtained  from  this 
reflection  study.  Over  the  years,  a  great  deal  of  effort 
has  been  expended  in  examining  crustal  properties  by  other 
means;  refracted  body  waves,  surface  waves,  potential  field 
methods  and  geochemical  observations.  With  so  much  interest 
focussed  on  the  earth’s  crust,  it  seems  appropriate  to  present 
a  brief  discussion  of  its  salient  features. 
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Usage  of  the  term  'crust'  to  describe  the  outer 
layers  of  the  earth  is  a  relic  from  ancient  beliefs  in  which 
it  was  thought  the  earth  was  a  thin  spherical  shell  enclosing 
a  hot  liquid.  However,  while  retaining  the  term,  modern 
usage  has  generally  redefined  it  as  being  that  part  of  the 

w/  /  ,  x 

earth  above  the  Mohorovicic  (M)  seismic  velocity  discontinuity 
(Steinhart  and  Meyer,  1961).  In  1909*  A.  Mohorovicic^  discovered 
this  seismic  discontinuity  from  an  analysis  of  earthquake 
waves  and  evidence  since  then  supports  its  universality,  giving 
depths  averaging  about  40  km.  under  continents  and  10  km. 
under  oceans.  Although  seismic  results  continually  verify 
its  existence,  the  nature  of  the  M  discontinuity  is  as  yet 
undetermined.  MacDonald  and  Ness  (i960)  have  proposed,  as  a 
consequence  of  phase  transition  theory,  that  it  may  be 
represented  by  a  gradual  transition  in  velocity  and  density 
from  basaltic  material  in  the  crust  to  eclogitic  material  in 
the  upper  mantle.  Their  work  indicates  that  this  transition 
zone  exists  over  a  distance  of  several  kilometers.  In  contrast 
to  this  viewpoint,  Nakamura  and  Howell  (1964)  have  suggested, 
as  a  result  of  analysis  of  refraction  arrivals,  that  the  M 
discontinuity  may  be  represented  by  a  relatively  sharp 
transition  extending  over  less  than  half  a  kilometer.  The 
observation  of  reflections  from  the  M  would  enable  the 
formulation  of  a  more  precise  theory,  be  it  either  of  the  above 
possibilities,  a  combination  of  the  two  or  the  proposal  of 
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something  quite  different.  Within  the  crust  itself,  there 
is  certainly  some  form  of  layering.  The  near-surface  sedi¬ 
mentary  system  of  layers  forms  a  relatively  thin  upper  region, 
but  until  1925  little  other  information  was  known.  In  that 
year  the  German  seismologist,  V.  Conrad,  discovered  evidence 
for  a  discontinuity  at  an  intermediate  depth  within  the  crust. 
Seismic  events  from  this  discontinuity  are  frequently  obtained 
in  continental  regions,  but  definite  proof  has  not  been 
obtained  that  it  is  a  world-wide  continental  phenomenon.  In 
general,  the  Conrad  discontinuity  is  interpreted  as  a  boundary 
between  an  upper  ’granitic’  zone  and  a  lower  1  basic’  zone 
within  the  crust,  the  nature  of  separation  being  unknown.  It 
should  also  be  pointed  out  that  on  seismic  records,  events 
from  the  Conrad  refractor  very  seldom  appear  as  first  arrivals 
so  that  identification  is  dependent  upon  the  quality  of  the 
data  and  subject  to  ambiguous  interpretation.  However,  if 
reflections  could  be  obtained  from  such  a  discontinuity  the 
ambiguity  would  be  removed  and  it  may  be  possible  to  determine 
whether  there  exists  a  widespread  continental  velocity  dis¬ 
continuity  within  the  crust.  Inherent  in  these  data  would  be 
information  about  the  nature  of  such  a  discontinuity. 

1.2  Seismic  Methods 

Mention  has  been  given  to  the  discovery  of  the  M 
discontinuity  by  A.  Mohorovicic  from  an  analysis  of  earthquake 
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waves.  The  recording  of  the  earth  motion  as  a  result  of  earth¬ 
quakes  and  the  subsequent  analysis  of  the  records  provided  the 
first  seismic  methods  for  ascertaining  properties  and  charac¬ 
teristics  of  the  earth.  For  years  this  was  the  only  seismic 
method,  but  the  unreliability  in  determining  the  origin  times, 
poor  knowledge  of  the  source  mechanism  plus  many  other 
problems  reduced  the  significance  of  the  results.  Steinhart 
and  Meyer  (1961),  pp.  3 point  out  the  shortcomings  inherent 
in  experimental  techniques  with  earthquake  seismograms.  With 
these  difficulties  the  natural  trend  was  toward  the  use  of 
controlled  explosions  for  studying  crustal  and  sub-crustal 
earth  properties.  Thus  were  born  refraction  and  reflection 
seismology,  although  the  initial  interest  for  crustal  studies 
lay  with  the  former  technique.  Long  range  refraction  surveys 
were  begun  in  the  1950’s  in  various  areas  and  contributed 
valuable  information  to  our  knowledge  of  the  crust.  However, 
because  the  energy  recorded  is  from  critical  refraction  along 
various  layers,  the  technique  cannot  identify  low  velocity 
layers  and  any  such  layers  contribute  to  errors  in  calculated 
apparent  vertical  velocities.  As  well,  the  calculation  of 
the  depth  to  any  refracting  layer  necessitates  an  assumption 
concerning  the  variation  of  velocity  with  depth  (Tuve  et  al, 
1954).  Some  of  these  difficulties  may  be  overcome  with  use 
of  the  seismic  reflection  method,  a  technique  whereby  energy 
propagating  downwards  in  a  near-vertical  direction  is  partially 


. 

,*v 

. 

*'.>*2  3  S  'lo  e  j.be.r  vcr  >i  iooq 

i  ••  b  r  '■  .  i< 

■ 

nl  nKj-9d  .dev; 

.  ?•  y  ::  1  ■ .  sv  d  ’t9ij?qq£ 

O  ‘j  o  ^  .  3 

9  IsrtSfW  p  U  9d  <1)0  + (10  i  DeX*)  9*1  0lflI3l98  9rfd  10 

xxy  sic  ,.ob  evn  1 


5 


reflected  from  velocity  discontinuities  of  sufficient  magni¬ 
tude.  In  addition  to  inherent  traits  of  greater  precision 
and  resolution,  the  method  allows  the  detection  of  discon¬ 
tinuous  low-velocity  layers  and  a  means  of  ascertaining  the 
average  velocity  to  any  reflecting  interface  (Green,  1938) • 

It  should  be  emphasized  at  this  point  that  the  ’sharpness'  of 
any  discontinuity  deduced  by  reflection  work  is  a  function 
of  the  wavelength  involved;  hence  the  unit  of  measurement  may 
vary.  The  reflection  technique  has  been  of  invaluable  use  in 
petroleum  exploration,  but  it  is  only  recently  that  attempts 
have  been  made  to  ascertain  characteristics  of  the  crust 
utilizing  this  method. 

1.3  The  Project 

It  has  been  noted  that  the  recording  of  near-vertical 
incidence  seismic  reflections  from  deep  within  the  crust 
would  be  of  tremendous  assistance  in  the  determination  of  the 
existence  and  of  the  nature  of  any  crustal  layering.  With 
this  thought  in  mind,  in  the  summer  of  1964  a  seismic  reflection 
program  was  begun  by  the  Geophysics  section  of  the  Department 
of  Physics  under  the  direction  of  Drs .  E.  R.  Kanasewich  and 
G.  L.  Cumming. 

The  field  operation  and  equipment  was  controlled 
by  certain  criteria  which  were  necessary  to  establish  the 
presence  of  primary  reflected  energy. 
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(1)  Since  such  energy  must  arrive  at  near-vertical 
incidence,,  the  direction  of  earth  movement  should  be  substan¬ 
tiated  with  the  use  of  both  vertical-  and  horizontal-component 
seismometers . 

(2)  As  previously  mentioned,  the  reflection  tech¬ 
nique  enables  the  calculation  of  average  velocities.  Thus 
enough  data  should  be  obtained  (preferably  about  a  common 
depth  point)  in  order  that  average  velocities  may  be  obtained 
from  reflecting  layers.  In  order  to  do  this,  the  quality  of 
the  data  must  be  sufficiently  adequate  that  correlation  from 
one  record  to  another  is  unambiguous.  Because  of  the  presence 
of  horizontally  travelling  energy,  this  generally  requires 
the  use  of  multiple  holes  and  geophone  patterns  (Section  2.4). 

(3)  Multiply-reflected  events  may  be  prevalent 
at  any  time  on  a  seismogram;  for  this  reason,  all  events 
should  be  regarded  as  multiples  unless  it  is  possible  to 
remove  any  ambiguity. 

(4)  Amplitudes  of  the  reflected  events  as  a  function 
of  distance  should  be  smooth  and  consistent  and  possess  some 
correlation  (or  valid  reasons  for  non-correlation)  with 
theoretical  amplitudes.  It  should  be  noted  that  near-surface 
velocity  irregularities,  shot  environment,  and  positioning 

of  the  seismometer  (on  soft  earth  or  hard  rock)  will  greatly 
influence  the  amplitude  on  actual  recordings. 
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(5)  Reflected  events  which  are  established  as  being 
from  some  velocity  discontinuity  within  the  crust  should  be 
obtained  at  different  recording  and  shooting  sites  within  some 
region,  providing  there  are  no  near-surface  noise  problems 
obliterating  all  events. 

(6)  Because  higher  frequency  wave  forms  are 
attenuated  more  rapidly  than  lower  frequency  ones,  it  would 

be  expected  that  the  frequency  content  of  any  deep  reflections 
would  be  considerably  lower  than  those  obtained  from  the 
sedimentary  layers.  This  criteria  thus  dictated  the  use  of 
low-frequency  geophones. 

The  experiment  was  begun  during  the  summer  of  1964 
in  a  region  in  southern  Alberta,  the  specific  location  chosen 
having  previously  been  used  during  the  recording  of  a  long 
reversed  refraction  profile  (Cumming  et  al,  1962).  By 
recording  a  reflection  profile  over  part  of  the  refraction 
profile,  it  was  hoped  that  the  reflected  events  could  be 
correlated  with  suggested  structure  from  refraction  work,  while 
at  the  same  time  yielding  more  precise  information  about  the 
nature  of  the  structure.  In  addition  to  a  continuation  of 
the  1964  profile,  work  in  the  summer  of  1965  included  the 
recording  of  profiles  at  two  different  locations.  Inherent 
in  the  entire  experiment  is  the  basic  and  corequisite  problem 
of  developing  a  technique  which  could  record  deep  crustal 


•  r-  '  b  ^3  :ooI:.v  emo* 

v  eelc  9t£t~  ;e-  ,  ci  on ci  >  ^nlbivoiq  tno.  3  n 

■  ■  .  u  -  .vo I  ino;  .  i:  00  •  . fo v, 

-?d  £n  •  .  IT  ■  ^fiemlJbsa 

,  inc  qc 

aor  no  '  ;0ol  0  llosqa  rid  tsv  *  IA  n'recifuoa  nl  no J,*3n  b  ni 
' 

0  t‘-  T-  '  .  r.  s  r  1  :;:'7  s  >3  tbnootn 

, 

^  ^  5  0;  ft-  '  o.  .  a'lijctc  J'ld?  arlcf  to  &iu3bci 

o&3  >  5  ■  .tqp If)  r'9b  to 


8 


reflections  continually  and  reliably.  It  is  felt  that  such 
a  technique ,  while  still  subject  to  much  improvement ,  has 
been  developed. 

In  conjunction  with  the  actual  field  data,  a  major 
part  of  the  author's  work  has  concerned  itself  with  the 
deduction,  construction  and  employment  of  analysis  techniques 
to  extract  useful  information  from  the  recorded  seismograms, 
or  provide  theoretical  data  which  may  be  compared  to  or 
applied  to  the  empirical  data.  To  this  end  computer  programs 
have  been  written  and  the  most  important  ones  have  been 
documented  in  the  appendices. 

Since  it  is  now  believed  the  reflection  technique 
will  provide  valuable  information,  future  plans  are  being 
made  for  the  continuation  of  the  program.  This  includes  an 
active  field  program  for  part  of  the  summer  of  19 66,  possibly 
extending  into  1967  for  completion.  More  time  and  attention, 
however,  will  be  directed  toward  gleaning  all  available 
information  from  present  and  future  data,  a  problem  which 
will  involve  extensive  data  analysis  techniques  and  theoretical 
model  studies.  By  means  of  projects  such  as  the  University’s 
reflection  program,  it  is  hoped  that  the  seismic  reflection 
method  will  become  a  proven  and  accepted  technique  for 
detailing  deep  crustal  structure. 
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1.4  Prior  Work 

Over  the  years ,  many  instances  of  deep  crustal 
reflections  have  been  reported,  but  many  of  these  were 
recorded  during  the  course  of  some  routine  seismic  investiga¬ 
tion  and  no  ancillary  equipment  or  techniques  were  used  to 
ascertain  that  they  were  primary  reflections.  The  first 
report  of  reflections  from  the  deep  crust  seems  to  be  that 
of  Junger  (1951) •  By  keeping  the  camera  running  to  10  seconds, 
reflections  were  observed  in  the  time  ranges  of  7-0  to  8.5 
seconds.  Junger  effectively  argues,  on  the  strength  of 
repeatable  results  and  identical  reflections  for  one  shot 
from  two  perpendicular  spreads,  that  the  energy  must  be 
arriving  from  some  reflecting  layers  inside  the  basement. 
Although  the  results  do  not  indicate  a  horizontally  continuous 
velocity  transition,  he  assumes  a  velocity  distribution  and 
suggests  the  range  of  depths  for  the  layers  to  be  l8  to  21  km. 
and  notes  that  these  depths  are  in  the  range  quoted  for  the 
granitic  layer. 

This  initial  publication  seemed  to  enliven  research 
and  subsequently  papers  were  published  in  many  countries. 

In  Germany,  Dohr  (1957)  was  one  of  the  first  to  publish  records 
including  late  arrivals.  His  work  is  monumental  primarily 
for  the  number  of  records  with  deep  reflection-type  events 
he  has  included.  With  so  much  data,  he  employed  a  statistical 
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approach  to  the  interpretation ,  showing  that  while  any  one 
reflecting  event  was  not  continuous,  on  a  statistical  basis 
there  existed  some  form  of  layering.  The  depths  correspond¬ 
ing  to  his  reflections,  assuming  some  velocity  distribution, 
range  from  a  shallow  layer  8  to  8.5  km.  deep  to  the 
deepest  at  30  to  31  km.  In  a  summary  on  crustal  structure 
in  western  Germany,  the  German  Research  Group  for  Exploration 
Seismology  (1964)  notes  other  areas  in  which  deep  reflections 
have  been  recorded  and  includes  a  bibliography  of  the  main 
researchers.  In  line  with  Dohr’s  work,  there  were  many 
records  available  for  analysis  and  a  statistical  approach 
was  again  the  basis  for  interpretation,  since  no  single 
reflection  event  could  be  correlated  continuously  over  more 
than  two  kilometers.  The  results  were  similar  to  those  of 
Dohr:  layering  existed  in  a  statistical  sense.  The  German 

researchers  argue  that  their  events  are  true  reflections, 
but  these  arguments  would  certainly  be  enhanced  had  they 
definitely  proven  the  events  were  due  to  primary  energy,  and 
that  the  wave-fronts  associated  with  the  reflection  arrived 
at  normal  incidence. 

Elsewhere  in  Europe,  B&th  and  Tryggvason  (1962) 
have  reported  on  the  first  seismic  investigations  of  deeper 
crustal  structure  in  Fennoscandia.  This  was  a  direct  attempt 
to  record  deep  near-vertical  reflections  but  the  results 
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were  certainly  not  conclusive.  The  authors  state  that  the 
arrivals  are  weak;,  of  erratic  occurrence  and  their  recogni¬ 
tion  is  based  on  phase  correlation*  while  admitting  that  some 
of  the  assumed  reflections  could  be  spurious  signals.  Never¬ 
theless*  they  expected  such  problems  and  only  identify  their 
events  as  'possible  reflections’.  Assuming  horizontal 
layers*  and  a  velocity  distribution  their  results  suggest 
a  depth  of  19  km.  to  the  Conrad  discontinuity  and  a  depth 
of  33  to  3^  km.  to  the  M. 

During  the  1950’s*  researchers  in  the  U.S.S.R, 
undertook  a  detailed  program  of  seismic  investigations  of 
the  earth’s  crust.  In  1960-1961*  Soviet  researchers  initiated 
an  experimental  program  to  record  deep  reflected  waves. 
Subsequently*  Beloussov  at  al  (1962)  reported  on  the  results. 
The  groups  involved  analysed  much  data  from  previous  seismic 
work  (primarily  supercritically  reflected  and  refracted 
waves)  and  adduced  an  optimum  recording  system.  This  involved 
the  grouping  of  receivers  and  shot  holes  to  enhance  desired 
signals;  figure  5  of  their  paper  emphasises  the  effectiveness 
of  their  system.  The  recordings  show  apparent  reflections  at 
times  as  late  as  18  secs,  and  greater.  In  addition  they  could 
phase  correlate  reflected  waves  along  their  120  km.  profile 
over  distances  varying  from  2-3  to  8-10  km.  Geophysicists 
in  the  Soviet  Union  have  recorded  many  deep  reflecting  events 
but  their  results  would  be  more  useful  if  they  indicated  those 
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events  which  could  be  due  to  multiply-reflected  energy,  as 
distinct  from  primary  energy.  However,  their  work  exemplifies 
the  amount  of  information  which  can  be  obtained  from  an 
organized  deep  reflection  program. 

In  North  America,  various  attempts  have  been  made 
to  record  near-vertical  incidence  reflections.  Subsequent 
to  Junger’s  report,  Shor  (1955)  presented  two  records  purported 
to  have  deep  reflections,  but  these  were  of  poor  quality  and 
the  possibility  of  multiply-reflected  events  could  not  be 
eliminated.  In  addition,  on  other  records  in  the  same  area 
no  reflections  could  be  found  because  of  high  reverberation 
levels.  Narans  et  al  (1961)  claim  to  have  obtained  reflections 
from  two  sites  in  northern  Utah.  Unfortunately,  their  records 
were  of  poor  quality,  requiring  a  statistical  approach  for 
interpretation,  and  they  could  not  exclude  the  possibility 
of  multiple  reflections.  Two  times  from  events  they  chose, 
when  combined  with  velocity  data  from  prior  refraction  work, 
gave  depths  corresponding  to  earlier  results  from  that  work. 
Robertson  (1963)  exhibited  a  record  section,  obtained  by 
extending  the  record  time  in  routine  petroleum  exploration  to 
about  six  seconds,  in  which  there  was  a  continuous  deep- 
reflecting  event  in  the  basement  complex.  The  reflections 
correspond  to  depths  much  less  than  those  reported  for  the 
Conrad  and  M  discontinuities  from  refraction  data.  Over  the 
29  km.  of  the  profile  in  southwestern  Alberta,  the  reflecting 
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layer  shows  considerable  dip,  increasing  in  depth  from  7.5  km. 
to  14  km.  Robertson  effectively  argues  that  the  reflections 
are  true  near-vertical  ones,,  but  he  only  speculates  on  possible 
interpretations.  More  recently  Dix  ( 1 965 )  has  reported  on 
reflection  studies  of  the  crust.  He  has  included  locations 
in  which  cross-records  were  shot  to  exclude  the  possibility 
of  interpretation  of  shallow  energy  from  a  side  reflector 
or  scatterer  as  primary  reflected  energy.  His  records 
exhibit  many  possible  reflecting  horizons,  with  a  strong 
sequence  of  pulses  centering  around  8.0  secs.,  corresponding 
to  a  depth  of  25  km.,  the  crustal  depth  for  the  region  of 
his  survey.  Finally,  Kanasewich  and  Gumming  ( 1 965)  have 
given  some  initial  results  from  the  1964  data  of  the  present 
project . 


'  t  r  )  ■  1 1  >  g  j  ct  c  . .  j* 


.  ■  '  ,  .  •  ;■  :  0- 


CHAPTER  2. 


EQUIPMENT,  LOCATION  AND  TECHNIQUES 

2,1  Instrumentation 

For  recording  purposes,  two  identical  sets  of 
equipment  were  arranged  in  two  panel  trucks.  This  equipment 
consisted  of  power  supplies,  low  frequency  amplifiers,  filter 
networks,  a  recording  oscillograph,  an  analog  magnetic  tape 
recorded  and  a  chronometer,  WWV  receiver  and  standard 
broadcast  receiver  for  timing  purposes.  In  figure  2.1  a 
block  diagram  shows  the  arrangement  of  the  apparatus.  A 
geophysical  contractor  was  hired  for  the  drilling  of  shot 
holes,  their  loading  and  detonation,  but  the  shot  instant 
was  recorded  in  a  suitably  equipped  university  vehicle, 

2,11  Amplifier  and  Recording  System 

Texas  Instruments  VLF-2  refraction  amplifiers 
were  used  throughout  the  operations.  Minor  modifications 
were  made  to  the  output  circuit  to  provide  the  necessary 
signal  amplitude  for  the  FM  analog  tape  recorder.  For 
rejection  of  60  cps  signals  induced  by  power  lines,  a  line 
balancing  switch  is  provided  with  each  amplifier.  In  addition 
to  this  switch,  useful  features  of  the  instrument  included  two 
input  plugs,  a  line  test  switch  to  check  that  proper  connec¬ 
tions  were  made  with  the  main  cable,  individual  step  attenua¬ 
tors  (6  db  steps),  and  a  24  db  per  octave,  6  position  high-cut 
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FIGURE  2.1  Arrangement  of  equipment  in  the  recording  truck. 
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filter  with  cutoff  frequencies  adjustable  between  8  and  48  cps . 
Recording  was  usually  done  with  a  filter  setting  of  24  cps. 

A  further  important  feature  of  the  equipment  is  dual  outputs -- 
high  gain  and  low  gain  with  a  deflection  ratio  of  4  to  1 . 

The  photographic  recorder  employed.,  an  RS-8U 
recording  oscillograph,  comprised  an  optical  system,  timing 
system,  recording  galvanometers,  a  paper  transport  system 
and  photographic  developing  box.  The  timing  lines  are 
obtained  from  a  100  cps  transistor  supply  and  a  synchronous 
motor  driving  a  drum  containing  slits.  For  recording  pur¬ 
poses,  thirty  galvanometers  are  available  of  which  twelve 
were  used  for  high  gain  traces,  twelve  for  low  gain  and  the 
remainder  were  used  for  timing  traces  or  ignored.  The  record¬ 
ing  was  made  on  roles  of  8"  x  200 1  seismic  photographic  paper. 
The  paper  speed  was  5"  per  second  allowing  a  timing  precision 
of  about  -  0.002  seconds. 

The  tape  recorders  utilized  were  Precision  Instru¬ 
ments  PS207  and  PS207A  instrumentation  tape  decks.  These 
allowed  six  channels  of  FM  seismic  information  plus  one 
channel  of  AM  time  information  (from  the  chronometer)  to 
be  recorded  on  l/2M  magnetic  analog  tape. 


I  ni  -■  w  1  1  O.;  fcnfi  tIbs  I  sir 

c  oX-9  )  oldqB'iaoloriq  br.o 


o  '3  ■  c  n:  le  .;:u  ,e  OOI  B  .oil  ten  r  i 


.  1  ^  ;  ■  \  to  ilj  t.  v Ini)  Tolom 

lo .'  ft  r  1.0  •  Cl..  V.  ■  3-19J  u  V  "7  1  ‘  ,  ?930q 


-  a? -ja  *1  sqBq  9fiT 


' ;  ■  "  ;• 


17 


2.12  Geophones 

Signals  from  six  of  the  recording  channels  were 
obtained  from  geophone  arrays  (Section  2.4)  consisting  of 
Electro-Tech  EVS-4  seismic  detectors.  The  main  characteris¬ 
tics  and  a  sensitivity  versus  frequency  plot  of  this  seismo¬ 
meter  is  given  in  figure  2.2a.  Also  used  to  record  earth 
motion  were  three  Hall-Sears  HS-10  detectors.  Figure  2.2b 
shows  the  principal  characteristics  of  these  geophones.  In 
some  instances,  signals  were  recorded  on  Texas  Instruments 
S-36  seismometers.  These  geophones  are  the  ones  employed 
for  the  University’s  refraction  experiments.  The  main  speci¬ 
fications  for  the  3-36  are  given  in  figure  2.2c. 

2.2  Location 

The  near-vertical  incidence  reflection  program  is 
being  carried  out  in  a  region  in  southern  Alberta  near  the 
town  of  Brooks,  The  main  reflection  profile,  called  the 
Lomond  profile,  is  located  on  a  twenty-mile  line  extending 
west  from  the  NEC  Sec.  24,  Twp.  17^  Rge .  18,  W4M  to  the 
NWC  Sec.  19,  Twp.  17,  Rge*  20,  W4M.  As  mentioned  before, 
this  location  was  chosen  because  it  had  previously  been  used 
for  part  of  a  reversed  refraction  profile  (Cumming  et  al,  1962) 
and  it  was  a  location  known  to  be  relatively  free  of  power 
lines.  All  the  reflection  data  from  the  summer  of  1964  plus 
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some  from  that  of  1965  were  obtained  along  this  line.  As 
well  in  1965*  further  data  were  acquired  from  a  short  profile 
just  west  of  Lake  Newell,  Alberta  (Newell  profile).  This 
profile  was  located  along  a  line  extending  west  from 
NEC  Sec.  22,  Twp.  17,  Rge .  15,  W4M  to  the  NEC  Sec.  22, 

Twp.  17,  Rge .  16,  W4M.  Records  Were  also  obtained  from  several 
locations  along  an  east-west  line  on  the  Defence  Research 
Board’s  Suf field  Experimental  Station.  Specifically, 
recordings  were  made  along  lines  extending  4  miles  east  from 
NEC  Sec.  33,  Twp.  17,  Rge.  9,  W4M;  4  miles  east  from  NEC 
Sec.  33,  Twp.  17,  Rge.  J ,  W4M  and  4  miles  east  from  NEC 
Sec.  31,  Twp.  17,  Rge.  5,  W4M.  Figure  2.3  shows  the  general 
area  in  wThich  reflection  work  is  being  done  and  the  positions 
of  the  profiles  including  the  designations  of  the  mile  positions. 

The  arrangement  of  holes  and  detectors  for  the 
Lomond  profile  was  an  adaptation  of  the  expanding  reflection 
spread  technique  (Musgrave,  1962)  with  center  point  at 
position  0  on  the  profile.  This  was  to  enable  average 
velocities  and  depths  to  be  calculated  for  any  reflecting 
horizon.  In  1964  the  centers  of  the  detector  arrays  were 
spaced  293  m.  (960  ft.)  apart.  Two  duplicate  sets  of  data 
were  obtained  as  for  each  position  of  the  detectors,  two  shots 
were  recorded,  there  being  a  separation  of  293  m.  between 
the  centers  of  the  hole  patterns.  The  duplication  of  data 
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FIGURE  2.3 

Figure  2.3  shows  the  seismic  reflection  profile 
locations.  Also  included  in  the  figure  are  some  well  loca¬ 
tions  from  which  data  were  obtained  to  compute  the  geological 
cross-section.  Listed  by  number  these  wells  are: 

2.  Cal -Standard  Parkland  4-12,  4-12-15-27  W4M,  1955 

3.  Grt.  Plains  Pageant  10-21,  10-21-18-21  W4M,  1962 

4.  Tenneco  Eyremore  10-15*  10-1 5-18-19  W4m,  1963 

5.  Cree  Scurry  Lomond  7-17*  7-17-l6-l8  W4m,  i960 

6.  Imperial -Cal st an  Lake  Newell  5-1*  5-1-17-W4M, 
I960 

7.  Cals tan  C.  P.  R.  Princess  No.  1,  13-22-20-12  W4M, 
1941 

8.  Richfield  Shell  Rapid  Narrows  11-20, 

11-20-16-4  W4M,  1957 
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was  to  facilitate  mixing  and  stacking  of  traces  at  a  later 
time.  Because  of  the  good  results  obtained  in  1964  and  in 
order  to  cover  more  distance  at  an  equivalent  cost,,  for  the 
1965  program  the  center-to-center  spacings  were  doubled  to 
586  m.,  except  for  some  special  shots. 

Only  two  shots  were  recorded  for  the  Newell  profile , 
these  being  detonated  at  positions  0  and  6  (figure  2.3)  with 
the  recording  trucks  located  at  positions  1  and  5.  A  spread 
of  1.609  km.  (l  mile)  to  either  side  of  each  vehicle  was 
utilized  with  a  spacing  of  5 86  m.  between  the  centers  of 
geophone  patterns.  No  attempt  was  made  to  employ  an  expanding 
spread  as  the  first  desire  was  to  record  reflections  in  a 
region  separate  from  the  Lomond  profile. 

Along  the  Suffield  line  a  number  of  shots  were 
made.  With  reference  to  figure  2.3*  patterned  shots  of  4 
holes  separated  by  30.5  m.  were  set  off  at  locations  0  and  4 
into  four  one-mile  spreads  between  these  positions.  Then 
the  spreads  were  moved  between  locations  12  and  16;  shots 
were  made  into  this  spread  with  identical  hole  patterns 
centered  at  positions  0*  4  and  16.  A  final  patterned  shot 
was  set  off  from  position  24  into  spreads  laid  out  between 
22  and  26.  For  all  profiles  it  was  necessary  to  lay  down 
about  25,000  feet  of  interconnecting  cable  and  install  212 
seismometers  for  each  position  of  the  recording  trucks. 
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2.3  Refraction  Summary  and  Geological  Cross-section 

The  Geophysics  laboratory,  University  of  Alberta 
has  previously  recorded  seismic  refraction  arrivals  across 
southern  Alberta  and  Saskatchewan  (Gumming  et  al,  1962  and 
Maureau,  1964) .  A  summary  of  the  data  obtained  from  these 
studies  is  illustrated  by  figure  2.4.  The  small  insert  shows 
the  location  of  the  profile. 

A  geological  cross-section  of  the  area  has  been 
computed  from  formation  tops  obtained  from  available  well 
log  data.  This  information  is  presented  in  figures  2.3  and 
2.5.  Digitised  acoustic  logs  from  the  Pageant  10-21  and 
Eyremore  10-15  wells  were  used  to  compute  the  upper  two 
average  vertical  velocities.  The  first  value,,  3*2  km. /sec.  , 
applies  from  the  surface  to  the  top  of  the  Mississippian; 
the  second,  5.8  km. /sec. ,  from  the  top  of  the  Mississippian 
to  the  top  of  the  Nisku  formation  in  the  upper  Devonian.  The 
two  basement  velocities  have  been  transferred  from  the 
refraction  summary. 

2.4  Geophone-hole  Patterns 

The  primary  information  desired  was  recorded  on 
arrays  of  EVS-4  seismometers  in  conjunction  with  hole 
patterns,,  the  combined  system  being  an  effective  filter  of 
long  period  surface  waves.  Six  channels  of  information  per 
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truck  were  obtained  with  this  system  and  the  signals  from 
these  geophone-hole  patterns  were  fed  into  the  FM  tape 
recorder  as  well  as  into  the  photographic  system. 

The  seismometer  array  utilized  was  a  ' tapered T 
array  of  l6  geophones.  This  consisted  of  twelve  seismometers 
being  laid  along  a  line  symmetric  about  the  take-out  location., 
with  a  spacing  of  12.2  m.  (40  ft.)  between  detectors*  result¬ 
ing  in  a  total  spread  length  of  134  m.  (440  ft.).  An  additional 
4  seismometers  with  the  same  spacing  were  laid  out  at  the 
center  of  the  system  to  give  the  tapered  array. 

The  normalized  response  of  such  a  system  of  geo¬ 
phones  (an  even  number  of  locations)  can  be  computed  using 
antenna  theory  (after  Savit  et  al*  1958) • 
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G  is  the  normalized  response 


a^  are  the  weighting  factors  for  2m  geophones 
located  at  distances  from  the  center  of  the 
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For  the  case  discussed,  then,  ai=a2=a3 
Dq  =  (12-1)40  =  440  feet.  Included  in 
of  (2.1)  for  the  designated  parameter 


=a^=a^=a^=l ; 

figure  2.6  is  a  plot 
values . 


The  hole  patterns  employed  were  quite  simple, 
principally  for  the  sake  of  economy.  In  1964  all  arrayed 
shots  except  the  first  two  were  made  with  a  pattern  of  4 
holes,  separated  by  6l  m.  (200  ft.)  and  symmetric  about  a 
central  axis,  each  hole  normally  being  loaded  with  50  lbs.  of 
explosive  charge.  The  initial  arrayed  shots  were  made  with 
separation  distances  of  30.5  m.  However,  in  1965^  an  attempt 
was  made  to  record  all  shots  with  a  hole  pattern  of  4  holes 
separated  by  30.5  m...  This  was  a  result  of  a  computer  program 
being  written  which  could  produce  many  response  curves, 
enabling  the  best  filtering  system  to  be  chosen. 


The  normalized  response  for  an  even  number  of 
holes  is  given  by  an  equation  identical  in  form  to  (2.1). 


H  = 


S»=l  bl  QOS 


(2.2) 


where 


H  is  the  normalized  response 

b.  are  the  weighting  factors  (relative  charge  size) 

v 

for  2n  holes  located  at  distances  from  the 
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o  o  e  y 


2,1 -1 


d' 


S  °  °  0  3 


±  2«zi  d, 


=  (2n-l)  dT  =  length  of  the  hole  pattern. 


Figure  2.6  Includes  the  hole  response  for  the  system  of 
4  holes ^  spaced  30.5  m.  apart.  The  total  response  of  the 
hole -geophone  combination  when  all  holes  and  geophones  are 
lying  along  the  same  line  is  given  by  the  products  of  the 
two  amplitudes.  This  curve  is  also  given  in  figure  2.6. 

Figure  2.7  gives  essentially  the  same  information  as  the 
previous  one  except  the  hole  response  is  for  4  holes  located 
6l  m.  apart.  In  addition,,  it  includes  the  response  curve 
of  an  infinite  number  of  geophones  on  a  line  over  the  distance 
of  134  m.  for  comparison  purposes.  It  will  be  noted  from 
figure  2.6  that  the  combined  response  is  an  effective  filter 
for  attenuation  of  arrivals  with  an  apparent  wavelength  less 
than  270  m.j  while  figure  2.7  indicates  an  effective  filter 
for  apparent  wavelengths  less  than  400  m.  Previous  mention 
has  been  given  for  the  choice  of  two  different  hole  pattern 
lengths.  Appendix  1  includes  a  derivation  of  the  responses 
for  odd  numbers  of  holes  or  geophones 3  and  some  additional 


response  curves. 
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FIGURE  2.6 


Amplitude  responses  of  geophone  and  hole 
patterns;  hole  spacing  is  30.5  m. 
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30 


patterns;  hole  spacing  is  6l  m. 
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As  previously  mentioned,,  use  was  made  of  Hall- 
Sears  geophones.  Three  of  these  per  truck  were  arranged 
at  one  location  with  three  take-out  leads  such  that  vertical , 
horizontal  radial  and  horizontal  transverse  components  of 
earth  motion  were  recorded  separately.  Component  recording 
of  the  ground  motion  for  deep  reflection  work  is  very 
important  for  the  output  of  these  geophones  indicates  the 
direction  from  which  the  energy  has  arrived.  Thus  while 
traces  from  the  geophone  arrays  may  indicate  the  arrival  at 
some  time  of  coherent  energy  over  the  spread  length 
(a  T ref lection T ) ,  the  directional  cluster  may  indicate 
whether  this  energy  has  arrived  from  a  truly  near-vertical 
di rection . 

In  general  the  third  type  of  geophone,  the  S-36, 
was  used  singly,  but  in  certain  instances  linear  arrays  of 
these  were  employed.  These  consisted  of  8  geophones  separated 
by  18.3  m.  (60  feet).  The  noise  attenuation  characteristics 
of  such  an  array  in  conjunction  with  the  hole  patterns  can 
be  calculated  in  a  manner  similar  to  the  previous  arrays. 
Appendix  1  includes  a  normalized  response  curve  characteristic 
of  this  array  and  hole  combination. 
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CHAPTER  3 


THE  DIGITIZED  RECORDS 

3.1  Method  of  Digitization 

Conversion  to  digital  form  of  the  analog  FM 
magnetic  tape,  recorded  along  the  Lomond  profile  in  1964,  was 
performed  by  the  technical  staff  in  our  geophysics  laboratory. 
Figure  3.1  is  a  block  diagram  of  the  procedure  followed. 

The  conversion  apparatus  was  originally  designed  for  the 
digitization  of  magnetotelluric  data  for  which  the  tape 
drive  speed  was  much  less.  This  necessitated  that  the  original 
tape  drive  speed  for  the  seismic  data  of  about  30  inches  per 
second  be  reduced  by  a  factor  of  sixteen  before  conversion. 

The  optimum  conversion  rate  for  synchronization  with  the 
paper  tape  punch  was  13.6  conversions  per  second,  machine 
time.  In  effect,  this  meant  217.6  conversions  per  second 
for  the  seismic  data,  yielding  a  digitizing  interval  in  real 
time  of  4.59  milliseconds. 

The  individual  channel  which  was  to  be  converted 
to  digital  form  was  selected  and  its  starting  time  controlled 
by  the  box  labelled  "clock  and  control".  The  filters  were 
specially  designed  to  prevent  aliasing  of  power  from  fre¬ 
quencies  greater  than  half  the  conversion  frequency.  Built 
into  the  filter  system  was  an  amplitude  modification  factor 
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Procedure  For  Analog  To  Digital  Conversion 
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which  could  be  set  at  +10  db..,  0  db. ,  or  -10  db.  depending 
on  the  dynamic  range  of  the  recorded  analog  signal.  Figure 
3.2  is  the  filter  response  plotted  as  a  function  of  real 
frequencies  as  opposed  to  the  actual  machine  conversion 
frequencies  which  are  a  factor  of  sixteen  less. 

AMF  =  Amplitude  modification  factor 
f  =  217.6,  the  data  conversion  freq. 

AMF 

AMF  -  20  db 

AMF  -  40  db 

AMF  -  80  db 


FIGURE  3.2.  Filter  response  as  a  function  of  real 


Conversion  frequency  in  real  time  (conversions 

per  second) 


conversion  frequency. 
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Signals  with  frequencies  greater  than  f  may  be  aliased 

c 

into  lower  frequencies  but  they  are  attenuated  by  40  db 
(a  factor  of  100)  or  more.  It  should  be  repeated  that 
prior  to  recording  on  tape,,  the  seismic  signals  above  24  cps 
were  attenuated  at  a  rate  of  24  db  per  octave. 

Once  all  channels  had  been  converted  to  digital 
form  on  paper  punch  tape,  the  data  were  transferred  to  punch 
cards  and  subsequently  put  onto  digital  magnetic  tape.  Since 
this  manipulation  of  data  could  be  a  possible  source  of  error, 
a  listing  was  carefully  scrutinized  to  ensure  that  all  data 
were  in  their  proper  sequence.  In  the  future,  such  a  tedious 
procedure  will  not  be  necessary  as  the  converted  data  will 
be  placed  directly  on  digital  magnetic  tape. 

3.2  Analysis  of  the  Digitized  Records 
3.21  Signal  enhancement  techniques. 

Subsequent  to  transferring  the  data  to  digital 
magnetic  tape,  a  computer  program  was  written  which  enabled 
the  data  to  be  plotted  on  the  IBM  1403  line  printer.  This 
was  a  quick  method  of  plotting,  enabling  a  visual  picture  of 
all  the  digitized  data  to  be  obtained  and  compared  to  the 
photographic  records.  Thus  a  correspondence  was  set  up  in 
which  each  digital  point  for  each  trace  could  be  related  to 
a  time,  accurate  to  within  5  milliseconds,  on  the  photographic 
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The  initial  step  in  the  analysis  of  the  9 6 
digitized  traces  was  the  normalization  of  amplitudes  and 
removal  of  the  normal  move -out.  For  the  former,  a  relatively 
simple  procedure  was  decided  upon:  one  representative  trace 
was  chosen  and  over  the  region  of  a  strong  reflection  at 
11.6  seconds  (Conrad  event),  the  amplitudes  of  all  other 
traces  were  normalized  to  that  of  the  representative  one. 

This  process  also  included  the  removal  of  any  D.C.  average 
prevalent  over  the  part  of  the  record  used  for  normalizing. 
Specifically,  a  half  second  interval  on  either  side  of  the 
Conrad  event  was  chosen;  the  mean  of  the  values  in  this  total 
interval  and  the  sum  of  the  absolute  value  of  the  deviations 
from  the  mean  was  computed  for  the  standard  trace;  the  same 
quantities'  were  computed  for  all  remaining  traces;  the  nor¬ 
malized  values  were  obtained  by  subtracting  the  mean  value 
of  each  trace  from  every  point  in  that  trace,  and  multiplying 
the  result  by  the  ratio  of  the  sum  of  the  deviations  for  the 
standard  trace  to  that  for  the  trace  under  consideration. 

It  should  also  be  pointed  out  that  while  the  total  record 
length  of  the  original  traces  varied,  the  output  of  the 
normalization  program  was  96  traces  of  equal  length.  Mathemati¬ 
cally  this  procedure  can  be  explained  as 
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2.  |  GStd  -  G8™ 
J  j 

2J  I  Gj  -«l 


(3.1) 


where 


and  the 


is  the  normalized  digital  value 

is  the  unnormalized  value 

j  refers  to  a  sum  over  the  1  second  interval 

T— T  refers  to  an  average  value  over  the  interval 

TstdT  superscript  refers  to  the  standard  trace. 


For  removal  of  normal  move -out ,  the  digital  value 
corresponding  to  the  bottom  of  the  ’picked’  trough  for  the 
Conrad  event  was  labelled  as  zero  point  on  each  trace. 
Effectively,  this  meant  that  the  times  of  events  representing 
reflections  from  the  Conrad  discontinuity  were  set  at  the 
reflection  time  (11.571  secs.)  of  the  representative  trace 
and  established  with  zero  moveout.  Any  events  at  different 
times  thus  had  the  move -out  of  the  Conrad  reflection  removed, 
so  compensation  for  the  normal  move -out  was  not  ideal.  In 
fact  a  calculation,  assuming  direct  ray  paths  and  utilizing 
known  refraction  data,  showed  that  the  difference  in  normal 
move-out  between  a  reflector  at  the  depth  of  the  Conrad  (34  km.) 
and  one  at  the  depth  of  the  M  discontinuity  (47  km.)  for  a 
horizontal  distance  from  the  shot  point  of  12  km.  was 
approximately  45  msecs.  This  horizontal  distance  corresponded 
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to  the  maximum  shot-to-detector  distance  for  the  1964  data. 

In  effect  then,,  for  record  times  greater  than  that  of  the 
Conrad  reflector  an  excess  amount  of  time  has  been  removed 
to  compensate  for  normal  move -out,  while  for  record  times 
less,  an  insufficient  amount  of  time  has  been  removed.  In 
the  future,  the  removal  of  normal  move -out  will  be  accom¬ 
plished  by  a  factor  which  is  a  function  of  recording  time  and 
horizontal  distance  from  the  shot  point.  Since  this  would 
involve  a  more  complex  computer  program,  it  was  deemed 
expedient  for  the  present  to  proceed  with  the  system  described. 
Figure  3.3  includes  two  traces  which  have  been  normalized, 
together  with  a  small  part  of  the  traces  to  the  same  scale 
before  normalization. 

As  a  check  on  the  zero  points  chosen,  neighboring 
traces  were  run  through  a  crosscorrelation  computer  program 
to  see  that  the  best  correlation  was  indeed  at  zero  lag. 

The  computer  program  is  similar  to  the  autocorrelation  part 
of  a  program  to  be  discussed  shortly.  Four  representative 
crosscorrelation  coefficients  are  plotted  in  figure  3*4  as  a 
function  of  point  lag.  In  conjunction  with  the  remaining  co¬ 
efficients,  they  indicate  that  the  proper  zero  points  were  chosen. 

It  has  been  previously  noted  that  reflections  for 
any  position  of  the  geophone  spreads  were  recorded  from  two 
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FIGURE  3-3.  Traces  representative  of  the  normalization 

and  stacking  of  digitized  data. 
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shot  points  separated  by  293  meters.  Hence  the  next  process 
was  to  add  equivalent  values  on  each  of  the  two  traces  which 
had  nearly  identical  travel  paths,  thus  reducing  the  total 
amount  of  data  by  one  half  to  48  traces.  The  result  of  this 
summation  process  on  two  traces  is  included  in  trace  5, 
figure  3.3. 


Stacking  or  compositing  of  traces  is  a  relatively 
direct  and  effective  means  of  increasing  signal -to-noise 
ratios.  Consequently,  a  computer  program  was  written  which 
would  stack  one  trace  with  its  neighbors  to  either  side, 
the  amplitudes  being  weighted  by  any  desired  factor.  Because 
of  changes  in  character  and  amount  of  step-out,  only  traces 
which  were  not  separated  by  a  large  distance  could  be  expected 
to  yield  a  significant  result. 

Two  major  runs  have  been  made  with  the  stacking 
program.  One  involved  the  successive  compositing  of  one 
trace  with  its  three  neighbors  to  either  side,  the  traces 
being  weighted  according  to  the  Gaussian  factors  (.0156, 

.0938,  .2344,  .3125,  .2344,  .0938,  .0156).  These  were  cal¬ 
culated  from  the  binomial  distribution  of  the  type 
B(j)  =  .  (i)  ^  ,  where  j  runs  from  0  to  6.  With  this  process 

the  original  48  traces  yielded  48  newly  filtered  traces. 

Trace  6  of  figure  3.3  was  obtained  in  the  manner  just  described. 
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The  second  major  run  involved  the  direct  compositing  of  six 
traces,  each  of  these  the  sum  of  two  traces  from  nearby  shot 
points,  to  yield  just  one  trace  (see  trace  7,  figure  3.3). 

With  this  procedure,  all  the  field  data  may  be  compressed 
into  eight  traces  representing,  in  general,  successively 
increasing  detector-to-shot  point  distances. 

Reproduced  in  figure  3-5  are  five  different  traces; 
the  first  three  are  a  result  of  600^  compositing  while  the 
final  two  are  a  product  of  the  Gaussian  stacking. 

It  will  be  noted  from  the  diagrams  that  the  stacking  procedure 
has  significantly  lowered  the  ’noise*  amplitude,  while 
retaining  the  original  strength  of  any  reflected  signals. 

Evident  on  these  records  is  a  lack  of  coherent  energy  before 
the  train  of  waves  beginning  near  11.5  seconds,  but  consider¬ 
able  energy  appears  to  be  present  later.  The  possibility  of 
multiple  reflections  is  clearly  indicated  in  figure  3*5  by 
the  duplicated  time  interval  between  prominent  events.  Un¬ 
fortunately,  the  stacking  procedure  has  not  notably  enhanced 
the  signal  amplitude  of  any  event  reflecting  from  the  M  dis¬ 
continuity.  However,  arrows  have  been  placed  in  the  figure 
around  15.0  secs,  depicting  coherent  energy  which  could  possibly 
be  reflected  from  M,  as  such  a  time  for  the  reflection  is  to  be 
expected  from  refraction  data.  That  the  pick  is  highly  ques¬ 
tionable  is  undoubtedly  true  and  future  work  both  in  the 
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recording  of  data  and  its  processing  will  include  as  one  of 
its  prime  objectives  the  identification  of  a  good  reflection 
from  the  M  discontinuity. 

3.22  Power  spectra  analyses 

Further  analysis  of  the  digitized  data  involved 
the  calculation  of  Fourier  transforms,  autocorrelation  coeffi¬ 
cients  and  autopower  spectra  for  selected  portions  of  the 
available  data.  The  necessary  equations  were  programmed  for 
the  University’s  IBM  7040  computer  and  all  calculations  were 
performed  on  it. 

The  Fourier  transform, 

+°° 

F(<“)  =  Ijp  /  f(t)  e~iaJt  dt  (3.2) 

-00 

was  resolved  into  real  and  imaginary  parts  with  the  integration 
being  based  on  a  linear  interpolation  between  digital  data 
points  (H.S.  Howe,  see  Appendix  2).  The  real  part  of  the 
transform,  for  programming  purposes,  is  then  given  by 

Re  [G(a>)]  =  -j-  2  [G^  )  -  G(t  )][cos  -  eos  ] 

CD  T 

(3.3) 

where  the  G(t^)  are  the  digitized  data  points  and  t  is  the 
digitizing  interval. 
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For  computing  Fourier  transforms,  an  interval 
representing  one  half  second  to  either  side  of  the  reflection 
at  about  11.6  seconds  was  chosen.  Seven  different  traces  from 
the  Gaussian  stacked  data  were  analysed,  these  representing 
varying  horizontal  distances  of  the  detectors  from  the  shot 
point.  Figure  3.6  is  a  plot  of  the  data  used  for  the  various 
computations  while  the  following  graphs,  figures  3.7  and  3.8 
give  the  moduli  of  the  Fourier  transforms  of  these  traces. 

Most  noticeable  on  these  figures  is  the  strong  amplitude  peak 
in  the  vicinity  of  10  cps .  Since  the  square  of  the  Fourier 
amplitude  is  representative  of  energy,  it  is  evident  that 
most  of  the  energy  in  the  traces  is  concentrated  within  a 
narrow  band  of  frequencies  extending  from  about  8  to  13  cps. 

This  is  presumably  due  to  the  strong  reflected  signa.1  while 
the  low  amplitude  side  peaks  are  probably  due  to  energy 
from  undesired  noise.  There  is  a  trend  for  the  moduli  to  peak 
at  higher  frequencies  with  increasing  distance  from  the  shot, 
a  point  which  will  be  mentioned  shortly. 

Power  spectral  calculations  were  made  following 
the  method  of  Blackman  and  Tukey  (1958),  (see  Appendix  2). 

The  digitizing  interval  of  t  =  .00459  secs,  gives  a  Nyquist 
frequency  (f^)  for  the  data  of  about  108  cps.  Let  N  be 
the  number  of  points  involved  in  a  computation,  L  be  a  correla¬ 
tion  lag  index  and  J  be  a  frequency  index  for  the  equally 
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spaced  spectral  calculations  in  the  range  0  <  f  <  f  .  The 
formula  for  computing  the  autocorrelation  coefficient  is 


C(L) 


1  yN-L 

N-L  i=l 


G.  T 
l+L 


(3.4) 


where  L  =  0,  1,  2,  .  .  .  ,  (N-l). 

For  calculation  of  the  autopower  spectra,  the  following 
formula  is  used: 

P(J)  =  2x  2^  W(L)  C(L)  cos  ^  (3-5) 

where  J  =  0,  1,  2,  .  .  .  ,  M 

M  is  the  specified  maximum  J  value. 

A  Daniell  window  (Appendix  2)  was  used  for  W(L),  a  weighting 
function . 

Autopower  spectral  calculations  were  made  for  the 
seven  traces  of  figure  3.6  and  for  a  one  second  interval 
preceding  the  data  shown.  In  all  cases,  the  number  of  points 
used  in  the  calculations  was  N  =  2l6,  the  constant  character¬ 
istic  of  the  Daniell  window  was  N/M  =  3  and  the  effective 
frequency  resolution  was  1.5  cps .  Figures  3. 9  to  3.12  are 
semi -logarithm  plots  of  the  frequency  versus  power.  No  correc¬ 
tion  has  been  made  for  instrument  response. 
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FIGURE  3.9 

Autopower  Spectra  for  I  Second  Interval  Preceding  Time  Interval  of  Figure  3.6. 
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The  graphs  of  the  first  two  figures ,  representing 
record  sections  with  no  evident  reflecting  event,  indicate 
considerable  energy  over  quite  a  wide  range  of  frequencies 
with  no  maxima  characteristic  of  one  frequency.  In  contrast, 
figures  3.11  and  3.12,  representing  record  sections  with  a 
highly  evident  reflecting  event  (figure  3.6),  show  graphs 
with  strong  maxima  at  around  11  cps .  The  energy  depicted 
by  such  peaks  must  thus  be  due  to  the  introduction  onto  the 
records  of  the  reflection  near  11.6  seconds.  Most  of  the 
energy  represented  by  the  peaks  is  concentrated  within  a 
range  of  frequencies  from  about  7  to  17  cps.  This  range  is 
slightly  greater  than  that  indicated  by  the  transforms  pri¬ 
marily  because  the  resolution  of  the  spectra  is  only  1.5  cps, 
so  that  the  energy  calculated  at  any  one  frequency  as  the 
Daniell  window  is  passed  over  the  autocorrelation  coefficient 
includes  energy  over  a  frequency  interval  of  1.5  cps.  As 
also  denoted  by  the  transform  plots,  the  peak  power  shifts 
from  about  10  cps  to  12  cps  as  the  shot-to-detector  distance 
increased.  It  is  believed  that  this  trend  is  due  mainly  to 
shot -hole  environment  and  the  surface  location  of  the  geophones 
The  correctness  of  the  supposition  should  be  resolved  with  the 
analysis  of  more  and  varied  data.  Because  of  stability,  and 
rejection  of  statistically  insignificant  changes  through  an 
averaging  process,  the  power  spectra  computations  are  to  be  pre 
ferred  over  the  simpler  Fourier  transform  calculations. 
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An  important  task  of  the  power  spectra  program  was 
to  give  an  indication  of  the  range  of  frequencies  over  which 
all  signals  are  received.  Thus  the  autopower  spectrum  was 
computed  for  traces  5  to  J  of  figure  3.3*  the  time  interval 
being  from  9*4  to  15.9  seconds.  With  N  =  1440  points  and  a 
characteristic  constant  for  the  Daniell  window  of  N/M  =  10, 
resulting  in  a  frequency  resolution  of  0.7  cps,  the  calcula¬ 
tions  were  made.  Figure  3.13  is  a  plot  of  the  resultant 
spectra. 

Since  the  graphs  represent  different  stacking 
techniques  including  the  same  trace,  a  few  comments  can  be 
made  concerning  the  results  of  stacking.  Both  methods 
involving  a  number  of  traces  reduced  the  amplitudes  on  either 
side  of  the  large  energy  peaks.  The  Gaussian  stack  reduced 
the  low  frequency  peak  most  and  increased  slightly  the  power 
of  the  center  peak  to  that  of  the  two  adjacent  peaks  in 
the  large  energy  region.  Reduction  of  the  amplitudes  at 
higher  frequencies  for  both  methods  was  about  the  same. 

The  most  important  point  from  the  spectra  is  the 
suggestion  that  band-pass  filtering  with  a  band  width  from 
5  to  15  cps  would  reduce  the  magnitude  of  the  side  bands 
which  are  thought  to  represent  undesirable  surface  waves  or 
wind  noise.  Following  this  suggestion,  future  plans  for 
analysis  include  digital  filtering  of  the  data  prior  to  any 
other  major  processing  techniques. 
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FIGURE  3.13  Autopower  Spectra  for  3  Traces} 

Time  Interval  of  9.5  to  15.5  Seconds 
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CHAPTER  4. 


SYNTHETIC  SEISMOGRAMS 


4.1  The  Synthetic  Seismogram  Computer  Program 

The  basis  for  the  computation  of  the  synthetic 
seismogram  is  a  means  of  calculating  reflection  responses 
for  a  system  of  plane  layers  in  which  the  velocity  is  a 
linear  function  of  depth  within  each  layer  (Berryman  et  al, 
1958) •  Assuming  plane  waves  at  normal  incidence  to  the 
layers  and  no  density  contrast,  Berryman  et  al  deduced  that 
the  reflection  coefficient  as  a  function  of  frequency, 
for  the  (j-l)th  layer  is 


where 


h-Yj-i  +  h-R  bjpjH+R^ 

J 

1  -R 

O-iH-i  -  (bj"bj-:P  +  H VitrO 


2  2 

[l-4tt>2/b  1] 

J 


-p.  ..In 


v  . 

J 


u-i 


(4.1) 


h  =  !VT5)/(zj+rzj} 


Vj+1  ,  are  the  velocity  and  depth,  respectively. 

at  the  end  of  the  jth  zone. 


Equation  4.1  was  obtained  by  solving  the  equation  of  motion 
in  elastic  solids  and  applying  the  conditions  of  continuity 
of  stress  and  displacement  at  the  boundary  of  each  layer.  The 
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or  last  layer  is  unlimited  and  its  reflection  coefficient, 
R^,  is  zero.  In  the  first  layer,  R^  represents  the  ratio 
between  the  reflected  wave  and  incident  wave  amplitudes  and 
consequently  is  the  transfer  function  of  the  entire  layer 
system  to  waves  at  normal  incidence. 

The  value  of  R^  is  computed  from  4.1  by  an  iterative 

procedure  beginning  with  R^  ^  (R^=0),  and  so  gives  the 

reflection  coefficient  for  plane  waves  of  normal  stress. 

For  displacement  the  calculated  value  is  multiplied  by  -1. 

Since  R.  is  indeterminate  if  either  b.  or  b.  is  zero, 
j-1  33-1 

the  equation  may  be  solved  for  this  special  case  or  more 
simply  for  computer  processing,  one  can  set  the  velocity 
gradient,  bj,  equal  to  a  very  small  value  such  as  0.001. 

The  latter  was  done  in  the  author’s  program. 

In  addition  to  having  the  reflection  response  of 
a  layered  system  as  a  function  of  frequency,  it  is  necessary 
to  have  some  form  of  an  input  to  the  system  before  a  synthetic 
seismogram  can  be  generated.  A  Ricker  wavelet  (Ricker,  1953) 
has  often  been  used  as  a  close  approximation  to  a  seismic 
pulse.  Examples  of  its  use  are  given  by  Jones  and  Morrison 
(1954)  and  by  Berryman  et  al .  The  form  of  the  Ricker  wavelet 
used  by  the  latter  consists  of  the  sum  of  the  fifth  and  sixth 
derivatives  of  the  error  function,  $(x)  : 
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G(x)  =  -  [b^(x)  +  <t>^(x)] 


(4.2) 


where 


(4.3) 


and  x  =  t/T,  t  representing  time  and  T  a  parameter 
characteristic  of  the  period  of  the  pulse.  The  synthetic 
seismogram  is  then  obtained  from  the  Fourier  transform: 


+oo  i<^>(t  +  z/v  ) 

S(t)  =  /  R  (w)  G(a))  e  °  d^  (4.4) 

-00 

where  G(^>)  is  the  complex  spectrum  of  the  incident  pulse 
and  z/vq  is  a  time  delay  factor.  The  Fourier  integral 
was  evaluated  by  the  method  described  in  Chapter  3  (or  see 
Appendix  2) . 


Since  it  was  imperative  that  an  accurate  working 
program  be  obtained.,  the  initial  development  of  the  program 
centered  about  a  duplication  of  the  results  of  Berryman  and 
his  group.  That  such  a  program  has  been  developed  is  evidenced 
by  figures  4.1  and  4.2.  The  former  is  a  copy  from  Berryman 
et  al  for  one  case  which  they  considered;  the  latter  is  the 
result  obtained  by  the  author.  In  similar  comparisons 
between  the  different  cases  which  Berryman  et  al  considered, 
the  author Ts  results  corresponded  nearly  identically,  small 
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FIGURE  4.1  After  Berryman  et  al  (1958). 
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differences  between  the  two  being  due  to  a  slight  difference 
in  the  period  of  the  input  pulse.  It  will  be  noted  that  the 
amplitude  scales  on  the  two  diagrams  are  different,  but  that 
corresponding  ratios,  such  as  peak-to-peak  reflection  to 
incident  amplitudes,  are  the  same. 

A  few  points  should  be  mentioned  concerning  the 
synthetic  seismogram  program.  The  time  involved  to  perform 
the  Fourier  integrations  is  very  long.  For  example,  a 
synthetic  seismogram  l6  seconds  long  took  7  hours  on  the 
IBM  7040  computing  system.  This  large  time  factor  was  mainly 
due  to  the  necessity  of  using  very  small  increments  in 
frequency  in  order  to  obtain  sufficient  resolution  to  delineate 
just  the  reflecting  events  and  not  obtain  spurious  events  of 
low  amplitude  caused  by  an  inadequate  superposition  of 
amplitudes  and  phases.  To  illustrate  this  point,  the  syn¬ 
thetic  seismogram  program  was  run  for  an  assumed  crustal 
model  using  a  frequency  increment  of  l/l6  cps .  The  resultant 
record  had  a  Conrad  reflector  but  an  extremely  weak  Mohorovicic 
reflector  plus  two  other  events,  at  about  6.4  and  10.1  seconds, 
which  could  not  be  explained  physically.  By  decreasing  the 
increment  to  l/40  cps  the  magnitude  of  both  the  Conrad  and 
Mohorovicic^  events  was  increased  significantly  and  no  un¬ 
explainable  events  were  found. 
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The  maximum  frequency  interval  which  can  be  used 
is  governed  by  the  length  of  record  section  desired.  For 
example,  with  an  interval  of  0.5  cps  the  maximum  record 
length  obtainable  is  2  seconds  because  the  entire  record 
for  the  two  seconds  repeats  continuously,  as  is  to  be  ex¬ 
pected  when  an  integral  form  is  reduced  to  a  series.  The 
inverse  transform  involves  the  terms  "cos  2tt  ft"  and  with 
f  =  0.5  cps,  this  term  has  the  same  value  for  t  =  0  as  for 
t  =  2.0  since  the  cosine  term  is  cyclic  with  a  period  of 
27 r  radians. 

Future  plans  include  a  modification  of  the  syn¬ 
thetic  seismogram  program  to  include  an  automatic  volume 
control  factor  and  possibly  some  function  representing  the 
instrument  response  which  is  included  in  a  true  field  record. 
Of  prime  concern  will  be  a  detailed  analysis  of  the  program 
to  see  if  its  operation  time  can  be  reduced,  as  this  is  a 
severe  limiting  factor. 

4.2  The  Mohorovicicf  Discontinuity  as  a  Transition  Zone 

One  of  the  first  points  investigated  was  the  range 
over  which  a  velocity  discontinuity,  represented  by  a  linear 
gradient  from  a  lower  to  a  higher  velocity,  could  prevail 
and  still  give  a  good  reflection.  Figure  4.3  includes  some 
results  from  this  study.  It  will  be  noted  from  the  diagram 
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that  a  change  in  velocity  from  7-32  to  8.25  km. /sec.  (repre¬ 
sentative  of  the  refraction  velocity  values  at  the  M  discon¬ 
tinuity)  over  a  distance  the  order  of  l/4  wavelength  gives 
a  distinct  reflected  pulse,  but  over  greater  distances  such 
a  distinctive  pulse  is  not  obtained.  The  width  of  such  a 
l/4  wavelength  transition  zone  would  be  l/5  kilometer  for  a 
frequency  of  10  cps  or  2  kilometers  for  a  frequency  of  1  cps . 
Therefore,  if  a  reflection  is  obtained  from  a  transition  zone, 
the  transition  must  be  2  kilometers  or  less  in  extent, 
dependent  upon  frequency  and  the  velocity  contrast.  A 
larger  reflection  coefficient  could  be  obtained  by  an  inter¬ 
ference  pattern  from  a  group  of  closely  spaced  transition  layers, 
but  the  general  principle  illustrated  above  still  applies. 

4.3  The  Synthetic  Seismogram  from  Acoustic  Log  Data 

For  the  purpose  of  making  a  synthetic  seismogram, 
digitized  acoustical  logs  from  two  wells  at  LSD  10-15-18-19  W4M 
and  LSD  10-21-18-21  W4M  (figure  2.3)  were  acquired.  The 
well  information  was  obtained  from  Petrex  Data  Agencies,  Ltd. 
on  IBM  punch  cards  and  consisted  of  the  travel  time  in  micro¬ 
seconds  over  a  one-foot  interval.  Since  velocity  and  depth 
values  are  required  for  the  synthetic  seismogram  program, 
this  necessitated  that  a  program  be  written  to  convert  the 
time  data  to  velocity  data.  This  was  done  but  within  the 
program  the  velocity  values  were  computed  as  an  average  over 
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a  ten-foot  interval.  Synthetic  seismograms  characterized  by 
input  pulses  with  dominant  frequencies  of  20,  28,  and  40  cps 
were  generated  from  the  well  log  data.  As  would  be  expected 
the  seismograms  with  the  higher  frequency  input  pulse  exhibited 
much  more  detail  than  the  corresponding  ones  with  a  lower 
frequency  pulse.  A  reasonable  correlation  between  the 
synthetic  seismograms  and  field  records  from  the  same  area, 
kindly  provided  by  Tenneco  Oil  and  Minerals,  Limited,  could 
be  made . 


The  principal  reason  for  generating  synthetic 
seismograms  was  to  obtain  some  indication  of  the  effects  of 
multiple  reflections  on  deep  reflecting  events.  Some 
examples  of  simple  multiple  reflections  which  could  be 
produced  are  illustrated  in  figure  4.4. 

In  making  the  synthetic  seismogram,  a  horizontally 
layered  model  of  the  crust  was  considered,  the  information  at 
depth  taken  from  the  refraction  data.  In  order  to  generate 
multiples  between  the  surface  and  any  strong  reflector  within 
the  sedimentary  strata,  a  near-surface  velocity  of  1.22  km. /sec. 
(4000  ft. /sec.)  was  assumed,  the  velocity  increasing  over  a 
3  meter  interval  to  2.90  km. /sec.  (9500  ft. /sec.),  this  system 
then  representing  a  velocity  step  function.  The  velocity- 
depth  values  from  the  acoustical  log  nearest  the  Lomond 
profile  formed  the  next  part  of  the  model,  these  representing 
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FIGURE  4.4  Examples  of  anticipated  simple  multiple  reflections. 
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the  sedimentary  section  to  a  depth  of  about  1.67  km.  From 
this  depth  to  34  km.,  the  model  included  a  slow  linear  increase 
in  velocity  from  about  5-90  to  6.30  km. /sec.  To  represent 
the  Conrad  discontinuity,  there  was  an  abrupt  increase  in 
velocity  at  34  km.  from  6.30  to  7-32  km. /sec.  over  an  interval 
of  3  meters.  The  M  discontinuity  was  represented  by  a 
sharp  increase  from  7-35  to  8.25  km. /sec.  over  3  meters  at  a 
depth  of  47  km.  All  of  the  induced  discontinuities  were 
included  essentially  as  a  velocity  step  function  to  obtain 
maximum  reflection  coefficients.  With  this  model,  nearly 
the  worst  case  of  multiple  reflections  may  be  studied. 

The  input  pulse  was  the  same  one  described  previously, 
but  with  a  dominant  frequency  of  13  cps .  In  figure  4.5* 
the  upper  half  of  the  diagram  shows  the  Ricker  wavelet  used 
as  the  input  together  with  its  Fourier  transform.  The  lower 
half  of  the  illustration  gives  the  reflection  coefficient  of 
the  crustal  model  as  a  function  of  frequency.  Within  the 
synthetic  seismogram  computer  program,  the  pulse  transform 
and  reflection  coefficient  are  multiplied  in  the  complex 
domain  and  the  inverse  transform  of  the  product  computed  to 
give  the  signal.  Plotted  in  figure  4.6  is  the  entire  16 
second  synthetic  seismogram  representing  reflected  energy  from 
the  crustal  model.  The  two  strong  deep  reflections  are  from 
the  Conrad  and  MohorovicicT  discontinuities .  Events  at  1.7  and 
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FIGURE  4.6 
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12.5  seconds  are  multiple  reflections  of  the  Interval  between 
the  initial  surface  event  and  the- strong  phase  at  0.9  seconds,, 
the  latter  representing  a  reflection  from  the  top  of  the 
Paleozoic  within  the  sedimentary  section. 

An  actual  field  recording  together  with  portions 
of  the  synthetic  are  presented  in  figure  4.7.  The  synthetic 
seismogram  indicates  that  multiply-reflected  events  are 
attenuated  quite  severely  within  about  1.5  seconds  after 
the  primary  reflection.  Comparing  the  field  record  with 
the  synthetic,  it  is  apparent  that  the  signal  on  the  field 
record  at  about  12.4  seconds  may  be  a  multiple.  However, 
the  event  at  13.8  seconds  is  probably  a  primary  reflection. 

On  the  basis  of  the  synthetic  seismograms,  it  has 
been  shown  that  strong  multiple  reflections  should  not  be 
a  problem  later  than  about  1.5  seconds  after  the  primary 
reflection.  This  is  an  important  conclusion  to  be  applied 
in  the  interpretation  of  any  field  recordings. 
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FIGURE  4.7*  Field  recording  and  selected  portions  of  the  synthetic  seismogram. 


CHAPTER  5. 


INTERPRETATION  OF  THE  ANALOG  RECORDS 

5.1  Suf field  and  Lake  Newell  Data. 

Experiments  were  carried  out  near  Lake  Newell 
and  on  the  Suffield  Experimental  Station  for  the  primary 
purpose  of  obtaining  deep  crustal  reflections  at  a  location 
different  from  the  initial  Lomond  site.  At  Suffield,  data 
were  obtained  along  three  6-kilometer  sections  of  an  east- 
west  line  nearly  33  km.  in  length.  The  signal -to-noise  ratio 
was  very  low  and  the  data  will  need  considerable  digital 
processing  before  it  can  be  interpreted.  The  noise  consisted 
of  surface  (Stoneley)  waves  generated  in  a  highly  weathered 
series  of  low  and  intermediate  velocity  lenses. 

A  portion  of  one  of  the  better  records  acquired  is 
presented  in  figure  5.1«  Coherent  events  which  possibly 
represent  reflected  energy  have  been  marked  on  traces 
recorded  through  use  of  arrays.  In  comparison,  the  single 
component  seismometers  show  a  considerable  amplitude  of  long 
period  surface  waves.  Although  no  attempt  has  been  made  to 
correlate  phases  across  the  various  Suffield  records  because 
of  their  poor  quality,  some  indication  of  what  layers  may  be 
reflecting  the  energy  can  be  given.  With  previous  results  as 
a  guide,  it  seems  likely  that  the  event  near  11.6  seconds 
represents  the  Conrad  reflector,  while  the  last  event  may  be 
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from  the  M  discontinuity.  Coherent  energy  at  times  prior  to 
11  seconds  could  be  indicative  of  the  low-velocity  layer 
which  is  suggested  in  the  refraction  summary  (figure  2.4). 

The  event  marked  at  13.8  seconds  is  found  on  a  number  of 
records  from  different  locations  and  seems  to  be  characteristic 
of  some  velocity-density  discontinuity,  not  mapped  by  refrac¬ 
tion  methods,  between  the  Conrad  and  Mohorovi^icf.  Only  when 
the  digital  processing  techniques  have  been  applied  will  it 
be  possible  to  give  a  definitive  meaning  to  some  or  all  of 
the  phases  marked  on  figure  5.1. 

The  small  program  planned  for  the  Lake  Newell  site 
was  severely  hampered  by  weather  conditions,  necessitating 
the  postponement  of  detonation  of  loaded  shot  holes  by  about 
2  l/2  months.  When  the  explosives  were  finally  discharged, 
not  all  holes  fired  and  consequently  the  effectiveness  of  the 
hole -geophone  patterns  as  a  filter  was  drastically  reduced 
(figure  5.2).  The  record  quality  for  the  third  set  of  traces 
for  which  all  shot  holes  discharged  is  generally  better  than 
the  others. 


The  shallowest  event,  near  10  seconds,  could  be  a 
reflection  from  the  low-velocity  layer  inferred  from  refrac¬ 
tion  data.  A  reflection  at  this  time  was  also  noted  at 
Suffield  and  on  a  record  from  the  eastern  end  of  the  Lomond 
profile.  The  event  at  11.6  seconds  is  probably  from  the 
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FIGURE  5.2„  Field  records  from  the  Lake  Newell  profile. 
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Conrad  discontinuity.  If  the  correlation  across  the  record 
section  is  correct,  then  there  must  be  approximately  l/2  km. 
of  west  dip  over  a  distance  of  9  kilometers.  Such  an  inter¬ 
pretation  agrees  with  gravity  data  which  indicate  a  high 
around  Lake  Newell  and  values  10  milligals  lower  at  Kitsim. 

5.2  The  Velocity  Determination  along  the  Lomond  Profile 

The  initial  phase  of  the  deep  crustal  reflection 
program  was  begun  along  the  line  now  referred  to  as  the 
’Lomond  profile’.  Detector  group  spacings  for  the  data 
acquired  in  1965  were  increased  from  the  293  meter  distance 
used  in  1964  to  5 86  meters.  This  proved  to  be  an  unfortunate 
decision  for  the  greater  detector  spacings  resulted  in  large 
step-outs  across  any  one  record  and  correlation  from  trace 
to  trace  became  difficult.  In  addition,  the  program  for 
1965  was  curtailed  by  several  weeks  of  rain. 

Approximately  half  of  the  data  recorded  along  the 
Lomond  profile  during  the  summer  of  1964  is  reproduced  in 
figures  5.3  and  5.4.  The  other  half  of  the  data  is  similar 
in  quality  and  represents  duplicate  coverage  with  a  slight 
offset  to  permit  application  of  signal  enhancement  programs. 
There  is  a  noticeable  lack  of  high  amplitude  coherent  energy 
from  4  to  11  seconds  after  the  shot  instant,  although  a  few 
weak  events  may  be  seen  on  some  traces.  The  first  prominent 
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FIGURE  5.3.  Lomond  profile  reflection  data  (1964)  between 

4  and  11  seconds  after  the  shot  instant. 
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FIGURE  p .  4 .  Lomond  profile  reflection  data  (1964)  between 

11  and  l6  seconds  after  the  shot  instant. 
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reflection*  near  11.6  seconds*  has  been  attributed  to  the 
Conrad  discontinuity.  The  arrow  indicates  the  particular 
cycle  which  was  chosen  as  the  zero  for  normalizing  the 
digitized  data  and  which  was  timed  to  compute  the  average 
vertical  velocity.  It  will  be  noted  that  although  no  event 
appears  to  exist  before  11  seconds*  there  is  a  suggestion  of 
coherent  energy  at  13.0  and  13.8  seconds  on  some  of  the 
traces*  the  latter  time  having  also  been  noted  on  other  records. 

Reproduced  in  figure  5.5  is  a  photograph  of  three 
different  field  recordings  for  the  time  interval  10  to  12.8 
seconds*  including  the  Conrad  reflector.  The  final  set  of 
traces  with  no  distances  marked  beside  it  are  the  low  gain 
equivalents  of  those  immediately  above.  One  point  this  set 
of  records  illustrates  is  the  widely  varying  signal-to-noise 
ratio  which  was  encountered.  The  upper  seismogram  is 
representative  of  a  region  where  the  signal-to-noise  ratio  is 
quite  low  while  the  middle  and  lower  recordings  are  character¬ 
istic  of  regions  with  high  ratios.  Traces  marked  fV*  designate 
a  single  vertical  component  geophone;  those  marked  TAT  designate 
arrays  of  seismometers.  A  comparison  of  the  two  sets  of 
traces  on  the  upper  record  gives  a  good  example  of  the  effec¬ 
tiveness  of  the  arrays  in  rejecting  long  period  noise.  To 
be  particularly  noted  on  the  lower  two  recordings  is  the 
definite  motion  of  the  vertical  component  seismometer  at  the 
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FIGURE-  5.5.  Lomond  profile  seismograms  (1964)  for  time  interval  10  to  12.8  seconds. 


82 


time  of  the  Conrad  event  while  at  the  same  time  the  horizontal 
component  geophones  (H)  show  no  unusual  activity.  This 
establishes  the  vertical  incidence  character  of  the  reflections. 

Prior  mention  has  been  given  to  the  fact  that  the 
program  of  1964  was  designed  to  utilize  an  expanding  spread 
profile *  this  having  approximately  a  common  depth  point. 

All  the  available  records  were  examined*  and  times  to  the 
Conrad  event  were  read  for  each  distance  from  the  center  of 
the  detector  group  to  the  center  of  the  hole  pattern.  The 
time  values  were  reduced  to  a  datum  plane  at  the  top  of  the 

Mississippian  system  at  an  approximate  depth  of  1.7  kilometers. 

2  2 

With  these  data*  an  X  versus  T  plot  was  computed.  From  the 
slope  of  the  line  and  its  y-axis  intercept*  determined  by  a 
least-squares  analysis*  the  average  vertical  velocity  from 
the  Mississippian  to  the  reflector  was  found  to  be  6.4  km. /sec. 
and  the  depth  of  the  reflector  was  calculated  to  be  35  km. 
below  the  surface.  The  standard  deviations*  assuming  no 
correlation  errors*  were  1.0$  in  both  the  velocity  and  depth. 
Comparison  with  the  refraction  data  of  figure  2.4  shows  good 
agreement  between  the  computed  values  and  those  indicated 
for  the  Conrad  refractor. 

That  the  correct  correlation  on  the  seismograms 
from  1964  has  been  made  is  evidenced  by  figure  5.6.  The  six 
stacked  traces  represent  J2  of  the  original  recorded  traces. 
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FIGURE  5.6.  IBM  1403  line  printer  output  of  six  &QG%  composited  traces 

The  circles  enclose  numbers  referring  to  the  trace  plotted 
(1  to  6).  The  output  also  includes  a  listing  of  the 
plotted  digital  values. 
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If  the  correct  phase  had  not  been  correlated,  the  similarity 
between  the  six  traces  would  not  be  as  striking  and  the 
pulses  would  not  be  so  well  defined.  This  illustration  also 
shows  the  form  of  the  output  from  the  plotting  program  which 
utilizes  the  IBM  1403  line  printer. 

A  continuation  of  the  expanding  spread  was  completed 
in  the  summer  of  1965  but  the  Bow  River  to  the  east  of 
position  +7  (see  figure  2.3)  limited  expansion  in  that 
direction.  Thus  the  additional  data  obtained  were  not  from  a 
common  depth  point.  Compounding  the  problems  already  mentioned, 
recordings  made  in  1965  were  most  generally  either  overshot 
or  undershot;  that  is,  the  gain  adjustments  of  the  amplifiers 
were  not  optimum  for  the  acquiring  of  good  visual  seismograms. 

Although  the  data  for  the  continuation  of  the 
Lomond  profile  were  poor,  it  was  found  possible  to  correlate 
the  reflection  associated  with  the  Conrad  discontinuity 
across  most  of  the  traces.  This  resulted  in  additional  time- 
distance  values,  which  were  added  to  those  from  1964,  and 
another  X  versus  T^  plot  was  computed  (figure  5.7) •  Prom 
these  data,  the  average  vertical  velocity  from  the  Mississippian 
to  the  reflector  at  a  calculated  depth  of  34  km.  was  found 
to  be  6.2  km. /sec.  Assuming  no  correlation  errors,  the 
standard  deviation  of  the  velocity  and  depth  values  was  0.4$, 
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However,  correlation  was  difficult  for  some  traces  and  this 

limits  the  reliability  of  the  results.  Ey  following  slightly 

different  phases  for  the  correlation  and  using  these  in  an 
2  2 

X  versus  T  plot,  it  was  found  that  the  depth  and  velocity 
magnitudes  varied  by  as  much  as  3 %  from  the  stated  quantities. 

As  well,  evidence  of  some  deep  structure  severely  restricts 
the  dependability  of  the  expanding  profile  in  this  region. 

5.3  Interpretation  of  the  Lomond  Seismograms. 

While  the  previous  sections  of  this  chapter  have 
shown  numerous  examples  of  reflections  from  the  deep  crust, 
these  have  only  been  noted  and  some  correlated  with  velocity 
discontinuities  known  from  previous  refraction  work.  In 
this  section,  the  results  will  be  consolidated  to  obtain  some 
structural  information.  A  Bouguer  gravity  survey  has  been 
made  and  this  tends  to  support  the  structural  implications 
of  the  reflected  events. 

Of  notable  significance  are  some  results  recorded 
during  the  1965  reflection  program.  These  consisted  of  several 
records  at  exactly  vertical  incidence  obtained  on  a  spread 
1320  feet  long  in  the  center  of  the  Lomond  profile.  A 
pattern  of  four  holes  with  only  10  pounds  of  explosive  per 
hole  was  detonated  and  all  records  showed  strong  deep  reflecting 
events.  One  of  these  is  illustrated  in  figure  5*8.  The 
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FIGURE  5.8.  Seismogram  at  vertical  incidence  from 

position  0  on  Lomond  profile. 
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strongest  event  is  undoubtedly  that  marked  just  after  12.0 
seconds,  although  weaker  arrivals  occur  at  11.5*  13.0,  13.8 
and  14 .8  seconds.  Utilizing  information  derived  from  the 
synthetic  seismogram,  the  13.0  second  event  may  be  a  multiple 
reflection,  but  it  is  not  likely  that  the  phase  at  13*8 
seconds  is  one.  Arrivals  at  the  latter  time  have  been  ; 
frequently  noted  on  other  seismograms.  The  last  event  at 
l4 .8  seconds  may  be  from  the  Mohorovicic  discontinuity. 

Figure  5.9  shows  the  strong  reflection  occurring 
at  12  seconds  on  the  vertical  incidence  seismogram  and  an 
equivalently  strong  reflection  at  11.6  seconds  on  two  nearby 
records.  It  seems  reasonable  to  correlate  these  two  events, 
but  to  what  is  the  time  discrepancy  of  nearly  one  half  second  to 
be  attributed?  It  is  possible  to  explain  these  results 
very  simply  by  assuming  a  high  angle  fault  deep  within  the 
crust.  Figure  5.10  illustrates  the  geometry  of  the  postulated 
fault  with  respect  to  the  seismic  ray  paths.  In  the  process 
of  faulting,  it  is  likely  that  any  horizon  in  the  vicinity  of 
the  fault  will  be  folded:  concave  up  on  the  deep  side,  concave 
down  on  the  shallow  side.  Because  the  amplitude  of  a  reflected 
wave  is  very  sensitive  to  curvature  of  the  reflecting  surface, 
a  focussing  effect  might  be  expected  to  enhance  an  event  from 
the  deep  side  of  the  reflector.  The  weak  event  at  11.5  seconds 
may  represent  diffracted  energy  from  the  edge  of  the  fault. 
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12  5  8  SURFACE 


FIGURE  5.10.  Geometry  of  the  postulated  fault 

with  respect  to  seismic  ray  paths. 
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Tracings  of  a  seismogram*  recorded  in  1965  about 
19  kilometers  from  the  explosion  source*  are  illustrated  in 
figure  5.11.  A  number  of  possible  reflections  exist  on 
this  recording.  That  many  of  these  come  from  near-vertical 
incidence  is  attested  to  by  the  strong  motion  of  the  vertical 
component  seismometer  and  the  weak  motion  of  the  horizontal 
component  ones.  The  phase  marked  C^*  at  about  12  seconds* 
has  been  correlated  with  the  Conrad  reflector  on  the  upthrown 
side  of  the  fault.  However  there  is  an  equally  strong*  if 
not  stronger*  event*  C^*  near  12.5  seconds*  a  time  which  for 
the  horizontal  distance  involved  would  be  consistent  with 
the  12  second  event*  representing  the  downthrown  side  of  the 
deep  horizon  on  the  vertical  incidence  recording.  The  long 
distance  seismogram  (figure  5.11)  also  contains  coherent 
energy  near  14.0*  14.5*  and  15.3  seconds.  Prom  the  synthetic 
seismogram  results*  it  is  not  likely. that  any  of  these  are 
multiple  reflections  of  earlier  events.  The  first  two  times 
could  be  related  to  reflections  from  the  shallow  and  deep 
sides  of  the  13.8  second  horizon  recorded  on  a  number  of 
records.  The  last  time  value  is  probably  indicative  of  an 
event  reflecting  from  the  shallow  side  of  the  M. discontinuity . 
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A  Bouguer  gravity  anomaly  map  for  southern  Alberta 
is  presented  in  figure  5.12.  On  this  map,  the  closely  spaced 
dots  represent  333  gravity  stations  at  which  readings  were 
made  by  the  Geophysics  group.  University  of  Alberta.  The 
remaining  dots  represent  stations  where  Bouguer  anomalies 
were  obtained  from  data  compiled  by  the  Dominion  Observatory. 
At  15  stations  where  both  university  and  government  data 
were  available,  it  was  found  that  Bouguer  values  from  the 
former  agreed  within  a  fraction  of  one  milligal  with  those 
listed  by  the  government. 

For  southern  Alberta,  the  regional  trend  is  negative 
towards  the  Rocky  Mountain  system,  but  this  trend  is  inter¬ 
rupted  by  a  pronounced  east -west  feature  centering  on  the 
Lomond  reflection  profile.  The  large  anomalies  are  likely  to 
be  due  to  density  variations  deep  within  the  crust.  The 
large  gradient  crossing  the  Lomond  profile,  marked  EAB, 
may  be  interpreted  simply  as  an  intrabasement  fault. 

Two  dimensional  model  calculations,  assuming  a 
throw  of  3  kilometers  down  to  the  north  and  some  plausible 
density  contrast,  have  been  carried  out.  These  gravity  com¬ 
putations  are  consistent  with  the  observed  gravity  field  if 
the  fault  penetrates  a  major  section  of  the  crust.  The 
postulated  fault  occurs  in  rock  of  the  Churchill  Geological 
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Province  (K-Ar  dates  of  1800  million  years)  but  does  not 
affect  Cambrian  rock.  Thus  the  suggested  east -west  fracture 
has  been  preserved  for  over  one  half  billion  years. 

More  information  will  be  necessary  to  verify  the 
structural  features  deduced  from  the  available  data.  The 
seismic  reflection  program  for  the  summer  of  1966  is  designed 
to  record  a  continuous  profile  across  the  postulated  fault. 

As  well,  further  gravity  data,  refracted  arrivals  from  the 
Precambrian  basement  and  ground  magnetic  observations  will 
be  obtained.  After  an  integration  of  this  geophysical  data, 
a  more  definitive  interpretation  of  deep  crustal  structure 
should  be  possible. 
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CONCLUSIONS 


Seismograms  recorded  in  southern  Alberta  clea.rly 
demonstrate  the  existence  of  near-vertical  incidence  deep 
reflecting  events  from  within  the  earth’s  crust.  These 
reflections  have  been  obtained  along  three  lines  over  a 
distance  of  150  kilometers.  Arrays  of  geophones  and  shot 
holes  have  been  designed  to  synthesise  a  filter  of  long 
period  surface  waves.  These  have  effectively  increased  the 
signal -to-noise  ratios.  Establishment  of  their  vertical 
incidence  character  has  been  achieved  through  the  use  of 
vertical-  and  horizontal-component  geophones. 

In  addition  to  photographic  seismograms,,  the  data 
were  recorded  on  PM  magnetic  tape  recorders  with  a  wide 
dynamic  range.  Digitization  of  these  analog  recordings,  their 
subsequent  normalization  and  the  compositing  of  appropriate 
traces  show  that  such  a  procedure  improves  signal -to-noise 
ratios.  The  application  of  more  sophisticated  data  processing 
techniques,  such  as  optimum  filtering,  deconvolution  and 
velocity  filtering,  before  stacking  should  enhance  the  ratios 
even  more. 


Other  analytic  techniques,  such  as  Fourier  trans¬ 
formations  and  autopower  spectral  analyses,  indicate  that  the 
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reflection  associated  with  the  Conrad  discontinuity  has  a 
dominant  frequency  of  about  11  cps.  Additional  autopower 
spectra  show  that  the  reflection  wavelets  have  the  greatest 
amount  of  seismic  energy  in  the  frequency  range  from  5  to 
15  cps.  Therefore,  additional  band  pass  filtering  in  the 
digital  domain  would  likely  improve  the  general  quality  of 
the  seismograms. 

For  theoretical  analyses,  a  synthetic  seismogram 
program  has  been  developed.  It  is  based  on  elastic  wave 
theory  in  a  layered  media  in  which  the  velocity  is  a  linear 
function  of  depth,  and  assumes  plane  waves  at  normal  incidence. 
Any  desired  seismic  wavelet,  in  the  time  or  frequency  domain, 
may  be  used  a3  the  signal  input.  Three  important  conclusions 
have  been  extracted  from  analyses  with  the  program: 

(1)  Multiple  reflections,  most  likely  between  the 
surface  and  the  top  of  the  Paleozoic,  may  be  prevalent  shortly 
after  any  major  reflecting  event  but  they  are  severely  atten¬ 
uated  within  1.5  seconds  of  the  start  of  the  primary  reflection. 

(2)  Any  transition  zone  from  within  the  deep  crust 
must  be  less  than  1  kilometer  in  extent  for  a  significant 


reflection  to  be  recorded. 
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(3)  Strong  reflections  on  the  synthetic  seismogram 
were  from  step-function  velocity  contrasts,  these  modelling 
the  deep  crustal  layers.  On  the  field  seismograms,  strong 
signals  reflecting  from  the  Conrad  discontinuity  have  been 
found  to  have  a  frequency  near  11  cps .  Together  with  the 
transition  zone  information,  these  facts  indicate  that  the 
extent  of  the  transition  zone  at  the  Conrad  must  be  small-- 
less  than  0.5  kilometer. 

A  major  advance  from  most  previous  data  is  the 
correlation  of  one  reflection,  attributed  to  the  Conrad  dis¬ 
continuity,  over  a  profile  30  kilometers  in  extent.  The 
expanding  spread  technique  enabled  the  average  vertical 
velocity  and  depth  to  the  Conrad  to  be  obtained  for  the 
first  time  and  formed  an  independent  check  on  assumptions  made 
in  refraction  analyses.  The  average  velocity  from  the  top 
of  the  Mississippian  to  the  discontinuity  at  a  calculated 
depth  of  3^  km.  +  3/  was  computed  to  be  6.2  km. /sec.  -  3 %. 

The  presence  of  deep  structure  beneath  the  region  has  limited 
the  precision  of  the  velocity  calculations. 

Reflecting  events  other  than  the  strong  Conrad 
reflector  have  been  indicated  on  a  number  of  different 
seismograms  from  different  locations.  One  of  these  is  an 
earlier  reflection  at  about  10  seconds  and  may  be  a  weak  event 
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from  a  low  velocity  layer  postulated  from  refraction  data. 

Two  later  reflections  around  times  of  13.0  and  13.8  seconds,, 
the  latter  having  larger  amplitudes,  frequently  appear. 

These  imply  that  additional  layering,  not  resolved  by  the 
refraction  method,  may  exist  within  the  deep  crust.  The 
latest  reflection  observable  on  any  of  the  seismograms  is  at 
a  time  near  15  seconds.  This  event  has  been  attributed  to 
the  Mohorovicicf  discontinuity,  its  arrival  times  agreeing 
consistently  with  times  to  be  expected  as  a  result  of  refrac¬ 
tion  information.  Definitively  proving  the  existence  of  any  of 
these  deep  horizons  requires  more  data  with  correlatable 
reflecting  events. 

Evidence  for  a  high  angle  fault,  with  a  throw  of 
3  kilometers  within  the  deep  crust,  has  been  presented.  A 
Bouguer  gravity  survey  indicates  a  large  anomaly  centering  on 
the  Lomond  profile.  Model  calculations  show  that  a  throw  of 
3  km.  is  sufficient  to  produce  the  necessary  gradient.  Dipping 
of  the  Conrad  reflector  to  the  west  as  suggested  by  seismograms 
near  Lake  Newell  is  also  accordant  with  the  observed  gradient. 

This  treatise  illustrates  the  wide  extent  of  informa¬ 
tion  which  the  reflection  method,  through  its  superior  resolu¬ 
tion,  can  add  to  knowledge  of  the  deep  crust.  It  is  to  be  hoped 
that  this  thesis  has  contributed  to  the  establishment  of  the 
seismic  reflection  technique  as  a  reliable  and  powerful  means 
for  deep  crustal  research. 
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APPENDIX  1. 

GE0PH0NE-H0LE  PATTERNS 


This  section  Is  an  extension  of  derivations  presen¬ 
ted  in  an  article  by  Savit  et  al  (1958)  entitled  "The  Moveout 
Filter".  Consider  an  array  of  N  seismometers  spread 
evenly  over  a  distance  D. 


- - D— - — — . . — 

-j|  d  V- 
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Since  the  form  of  any  wave  impinging  upon  the  seismometers 
will  be  a  function  of  both  time  and  distance,  the  input  to 
the  geophone  group  can  be  expressed  by  a  double  Fourier  series. 
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J=1  P^l 


[sin 


(oj  t 

v  P 


+ 


2tt 

h 


x  + 


0,  )] 
3,  P 


where 


=  angular  frequency 
A  .  =  wavelength 

J  <■>  A 

v  .  -  ■  =  apparent  phase  velocity  of  any  Fourier 


component . 


I. 


For  the  following  analyses,  a  single  frequency  and  wavelength 
is  considered  as  the  input  to  the  geophone  array. 

Ojr 

yin  =  A  sln  (a)t  +  —  x  +  <0  (Al.l) 

The  normalized  output  from  the  array  is  the  sum  of  individual 
outputs : 


y out  =  y  =  ^  A  sin  +  l2Lnd  + 


N  n-0 


A 


D 


It  is  evident  from  the  diagram  that  d  - 


1  2N-1  .  ±  fa)t  ,  ££ iii  u  QS 

y  N  n=0  v  N-l  X  J 


2mr  D 


(A1.2) 


We  utilize  Lagrange’s  trigonometrical  identity: 
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2m=0  Sln  +  my^ 
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and  A1.2  becomes 


y  = 


/Ntt  D\ 

sin 

N  sin  (—■  £) 


sin  (^t  +  +  $) 


(A1.3) 
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Comparing  A1.2  with  A1.3  we  note  that  the  output  of  the 
group  Is  equivalent  to  that  of  one  seismometer  placed  at 
x  -  D/2  and  having  the  amplitude 


.  /  Ntt  D\ 

Sin  ^N-l  A^ 

N  sin  (  7r 


D 

N-l  A 


We  wish  to  evaluate  the  effect  of  the  grouping  for 
noise  suppression.  Evidently  the  most  effective  would  be 
an  infinite  number  of  seismometers  spaced  over  the  distance  D. 
Thus  we  consider  the  limit  as  N  approaches  infinity  of  A^. 


lim 


sin 

~ttT) 

A 


sa 


00 


(A1.4) 


This  formula  was  included  in  the  computer  program  to  give 
the  equivalent  response  of  an  infinite  number  of  geophones 
over  a  distance  D. 


Consider  a  symmetric  array  of  2n  +  1  detectors  of 

ID  2D 

arbitrary  sensitivity  placed  at  locations  x  -  0,  ±—  br.»  7p 
.  .  . ±75-  with  sensitivities  given  by  a^,  a^ ,  aQ„  .  .  .  „  an» 

We  specify  that  the  output  of  the  array  be  unity  for  a  signal 
arriving  normal  to  the  surface,,  corresponding  to  a  signal 


. 
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with  infinite  phase  velocity  and  thus  infinite  wavelength. 
The  input  signal  is 

Ptt 

yin  =  cos  x  +  $)  (A1.5) 

The  output  of  the  group  becomes  with  d  =  ^  §• 


2tt 


y  =  Z^  +’n  a,  cos  (cut  +  kd  +  P) 
Jout  k=-n  k  v  A  ' 


(ai.6) 


Expand  the  cosine  term  within  the  summation  sign  of  A1.6. 


Ptt 

cos  (<^t  +  —  kd  +  $)  = 


O  'TT*  O  rTT~ 

=  cos  (u>t  +  <l>)  cos^ —  kd  -  sin  (,o>t  +  P)  sin  kd 


A 


However*  the  second  factor  will  not  contribute  to  the  summation 
since  the  a^  are  symmetric  about  x  =  0  and  the  sine  function 
is  odd.  That  is 

2iT+n  av  sln  f^kd  =  0 

k=  -n  k  A 

Then  using  the  symmetry  of  a^.*  A1.6  becomes 


'  .  I  '  '  i  t  .  O  .  . 


' 


sonls 


■r;  e.-  4X1  sxj  n  -.r(T 


107 


yout  =  ao  cos(a)t+ct>)  +  2  ak  cos(a)t+ct))  cos^ykd 

=  cos(u>t+4>)  [a  4-  2  2^__1  cos— y-  kd] 

=  G^  cos  (^t  +  b) 


wh 


ere  G-^,  the  amplitude  of  the  response  is 


~  .  n  2tt  n  , 

G-,  =  a  +  2  E.  _  a,  cos  — =r-  kd 
1  o  k=l  k  A 


Normalized  to  the  desired  unit  response  gives  the  amplitude: 


G 


a0  +  2  zk=i  ak  cos  vr  kd 

ao  +  2  2k=l  ak 


(A1.7) 


If  we  now  let  N  =  2n  +  1,  the  number  of  detectors,  then 


N-l  .  2  D  D 

n  =  -7s~  so  d  - 


N-l  2  N-l  # 


With  these  substitutions,  A1.-7  becomes 


G  = 
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(A1.8) 
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A  similar  analysis  will  give  for  M  =  2m  seismometers  located 
at  ±75-  ,  ±  rr’,  . ..,  ± ^  with  sensitivities  of  a^,  a^, 

.  ..,  am,  an  amplitude  response  of 


G  = 


M 

2k=l  ak  cos 


7 T 


M 

2,2  n  a, 
k=l  k 


2k-l  D 
M-l  A 


(A1.9) 


Equation  A1.9  is  identical  to  equation  2.1.  The  effects  of 
a  hole  pattern  can  be  derived  in  an  identical  manner.  If 
we  call  the  amplitude  response  of  the  hole  pattern  H,  then 
the  combined  effect  of  the  geophone-hole  pattern  is 


R  =  |  G  |  |  H  |  . 


(A1.10) 


A  computer  program  was  written  which  evaluated 
G,  H  and  R  for  any  desired  symmetric  configuration,  at  the 
same  time  calculating  the  response  of  an  infinite  number  of 
seismometers  over  the  specified  distance  D.  It  might  be 
noted  that  the  argument  of  tbe  cosine  term  for  the  geophone 
pattern  and  the  hole  pattern,  —  and  —  ,  respectively, 
will  in  general  be  different.  The  evaluation  of  the  combined 
response  must  be  done  for  various  wavelengths,  but  the 
individual  terms  of  G  and  H  must  be  considered  at  the 
same  wavelength.  For  this  reason,  two  parameters 
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„  DG  (NG  -  1}  dG 
G  A  “  X 

„  DH  %  -  ^  dH 
n  ~  A  ~  A 


were  used  and  Zu  was  derived  as  a  function  of  Zn,  specifying 
A  to  be  the  same.  This  gave 


(Nh  -  1)  dH 
(NG  -  n  dG 


(Al.ll) 


where  dTT  and  d~  are  the  distances  between  individual 

n  u 

holes  and  detectors,  respectively;  and  are  the  number 

of  holes  and  detectors,  respectively. 

The  following  graphs  show  the  amplitude  responses 
for  the  various  configurations  indicated  thereon. 
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FIGURE  A1.2 
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FIGURE  A1 . 3 
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FIGURE  A1.4 
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APPENDIX  2. 


FOURIER  TRANSFORMATION  AND  POWER  SPECTRA  ANALYSES 


A2.1  Fourier  Integration  on  a  computer 

The  Fourier  and  inverse  Fourier  transforms  are 

given  by 


S(^)  =  hr  S°°  e  ia)t  dt 

—00 

(A2.1) 

00  * 

s(t )  =  /  S(cu)  eia)t  d<u 

(A2.2) 

-00 


In  this  section,  the  equation  suitable  for  application  on  a 
computer  will  be  derived  for  A2.2.  The  method  is  analogous 
to  that  described  by  Howe  (1956)  for  equations  of  the  type 
A2.1. 


The  integral,  A2.2,  is  divided  into 
imaginary  parts. 


real  and 


00 

s(t)  =  /  [S^fco)  +  i  ST(^)][cos  <^>t  +  i  sin  ^t]  d^ 

-oo  K  1 

00 

=  /  [S  (u>)  cos  cut  -  ST  (cu)  sin  cot]  d^ 

-oo  R  1 

00 

+  i  /  [Sp(co)  sin  cot  +  Sj’(co)  cos  cot]  dco 
-00 


(A2.3) 
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Empirically,  the  signal  s(t)  must  be  real  so  the  imaginary 
part  of  s(t)  must  vanish.  We  note  this  will  occur  when 
SR(^)  is  even  and  Sj(^)  is  odd.  With  these  conditions 
A2.3  becomes 

00  00 

(A2.4)  s(t)  =  2  f  S  W  cos  0>t  d03  -  2  /  S]-(^)  sin  o>t  d oj 

-00  -00 

Empirically,  Snt^)  and  ST (^)  are  known  at  discrete  intervals 
of  frequency.  We  approximate  a  continuous  function  by 
assuming  each  varies  linearly  between  interval  points. 

Thus  we  first  wish  to  evaluate  each  of  the  above  integrals 
when  they  are  integrated  between  finite  limits,  assuming 
linear  variations  between  limits. 

Consider  first 


(A2.5) 


/“  S0(to)  cos  ait  dcu 
J  0  Rv  ' 


Let  and  be  the  values  at  the  beginning  and  end  of 

j  J+l 

the  (j+l)th  interval,  respectively.  The  corresponding  function 
values  are  and  ST3(cb.in).  Then  the  contribution  to 

n  j  K  j+l 

A2.5  from  this  interval,  assuming  the  linear  variation,  is 
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(A2.6) 


Now 
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cDt  do  =  i  [sin  (o  t)  -  sin  (o  t)]  (A2.7) 
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Using  A2 . 7  and  A2.8  in  A2.6  leads  to 


SrK+1)  -  s  (0)  ) 

ASr  =  - - - - - p— [cos(u>. 4lt)  -  cos(^.t)] 

R  (Vl-^)t2  J+1 

+  ^  SR^J+l)CUj+l  Sin(Vlt}  -  SR(Vl)CUJ  Bln(V} 


-  VWi  sin(Vit)  +  sr(cV“j  sin(Ut} 
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Spf^ui  )“,•  sin  (u> .  ,  t )  +  S„(<bu1  )a),  sin(co.t) 


R  J+l'  j 


(V  -  V} 


J+l 
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Simplifying  we  obtain 


S  fa)  )  .  s  (cd  ) 

ASr  =  — — ^ [cos(u>  -,t)  -  cos(^.t)] 

+  sR(y)  sln^J+it)  -  W  Bln(<aJt) 

t 


(A2.9) 


Summing  consecutive  intervals  to  obtain  the  total 
value  of  the  integral  A2.5,  we  notice  that  [S^f^ . . n ) sin (w . , _ t) ]/t 

K  J+l  J'r  l 

for  the  first  interval  equals  [SR(^. ) sin (^.t )  ]/t  for  the 
second  interval;  hence  these  terms  cancel.  If  we  sum  for  all 
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intervals,  all  interior  points  vanish  leaving  the  first  and 

last  points.  However,  at  03  =  0+  and  at  some  higher  frequency 

(depending  on  the  maximum  frequency  content  of  the  signal) 

the  values  of  SR  approach  zero.  That  is,  S^(o3  )  is  zero 

for  the  first  interval;  S^(o)  )  is  zero  for  the  last 

n  J+l 


interval . 

Thus  the  second  term  of  A2.9  vanishes  upon 

summation 

over  all  intervals. 

This  gives  for  A2.5 

/qSr(co)  cos  tut  dtu 


(A2.10) 

Sd(o).l1  )  -  SD(o).) 

=  2.  — — ^±1 - § — i  [cos  (03  t)  COS  (03. t)l 

J  (Vi-  “9  * 

where  the 

sum  is  over  all  the  intervals. 

identical 

Evidently  /q  S^(o3)  sin  o>t  d<^  will  give  a  result 

in  form 

/“  S_  (to)  sin  tut  dtu 


(A2.ll) 

ST(to  )  -  ST  (to.) 

=  2.  — — itl - A— ^  [sin  (to  t)  -  sin  (to  t)] 

3  (Vi-  *  3+1  3 

Using  A2.10  and  A2.ll  in  A2.4  allows  s(t)  to  be 

calculated.  The  form  of  the  equation  is  appropriate  for  use 
on  a  digital  computer. 


tBlsvxstfn'i 
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A2.2  Power  Spectra  Analyses 

The  theory  presented  In  this  appendix  is  a  modifi¬ 
cation  of  information  presented  by  Lee  (i960)  and  Blackman 
and  Tukey  (1958) . 

Assume  we  have  recorded  a  function  f(t)  generated 
by  a  random  stationary  Gaussian  process.  The  integral  of 
the  square  of  f(t)  must  be  finite.  As  well,  the  signal 
must  have  no  D.  C.  value: 


(A2  ..12) 


f 


lim 


T-*°° 


T/2 

/  , 

-T/2 


f(t)  dt  =  0 


The  variance  is  defined  as 


,  T/2  _  2 

C(0)  =  lim  4  I  U(t)  -  U  dt 
1  1  -T/2 


or  by  A2.12 


(A2.13) 


'2 


T/^ 

C (0)  =  lim^  i  /  .  [f(t)]  dt  . 


-T/2 


The  autocovariance  of  the  function  f(t)  is  defined  by 


(A2.14) 


T/2 

C(t)  =  lim^  /  .  f(t)  f*(t  +  t)  dt 


-T/2 


where  the  refers  to  the  complex  conjugate  and  t  is  the  lag. 


S'  ♦  ‘  A 
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In  this  study,  we  will  consider  real  signals. 

The  autocovariance  is  often  called  the  autocorrela¬ 
tion  coefficient  although  strictly  speaking  the  latter  term 
should  refer  to  the  normalized  function 


$  (  t  )  =  C  (  t  )  /  C  ( 0 ) 

The  signal  f(t)  in  the  autocorrelation  coefficient 
extends  over  the  range  -T/2  to  +T/2  but  f(t  +  t)  extends 
from  (-T/2  +  t)  to  (t/2  +  t).  To  avoid  complications, 
we  require  that  f (t)  =  0  if  t  exceeds  the  range  |  t  |  >  T/2. 

The  Fourier  transform  of  f(t)  is 

-j-00 

(A2.15)  F(“)  =  4=  /  f(t)  e'1(Ut  dt 

-00 

while  the  inverse  transform  is  given  by 

-{-00 

(A£.l 6)  f(t)  =  /  F(o>)  elQ)t  do) 

—00 

Substituting  for  f (t+  t)  in  A2.14  gives  for  the  autocorrela¬ 
tion  coefficient 

i  /T//2  f(t)  /MF(cu)eiC0(t  +  t) 

T-“  T  -t/2  -co 


C(t)  =  lim, 


d^>  dt 


•>  o  -  £  $&n.:  oo o  *r  f 
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or  changing  the  order  of  integration 

H"1  /  P  4-oo 

C(t)  =  lim^  i  /  F(co)  e1CUT  f(t)  ela)tdt  d<u 

Prom  A2. 15  with  <^>  =  -d) 


(A2.17) 


2ir  T/2 

C(t)  =  limT^,  /  .  F(co) 


-T/2 


iOf  ,  . 

e  F(-O)) 


We  now  define 

P(“)  =  limT^  f1  tF(co)  •  F(-co)] 
or 

Ptt  2 

(A2.18)  P(<»)  =  lim^  —  |F(cu)  I 

Equation  A2.17  becomes 

-|~°0 

(A2.19)  C(t)  =  /  P(o>)  elLjT  do) 

-00 

P(d))  is  the  power  spectrum  and  is  given  by  the  inverse 
Fourier  transform 

+°° 

(A2.20)  P(<0)  »  |=  /  :  ;  C ( T )  e‘iaJxdT 

-00 

Since  the  autocorrelation  function  is  even  this  can  be 
rewritten  as  a  one-sided  Fourier  cosine  transform 


*:  '  j-jsrtasjn:  U  i  to  ertt  snJtamrrio  to 


T  «kT 


•}  .  fpi.  ' 
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(A2.21)  P(^)  =  —  /  C(t)  cos  dx 

0 


It  should  be  noted  that  the  autocorrelation  function  and  the 
power  spectrum  of  a  random  function  are  related  to  each  other 
by  a  Fourier  cosine  transformation. 


In  reality.,  the  records  we  examine  are  finite  in 
extent.  For  such  cases  we  cannot  estimate  C(x),  the  auto¬ 
correlation  function  for  arbitrarily  long  lags.  Thus  we 
define  an  apparent  autocorrelation  function 

r  i(  t  - 1  t  |) 

.  f 

-|(Tr  -  I  T  I) 


(A2.22) 


Ca^ 


T  - 


X 


(t  -  -|r)f(t+  i-x)  dt 


with  I  t  I  <  t  <  T  where  t  is  the  maximum  desired  lag 

value  and  T^  is  the  record  length.  The  autocorrelation 
function,  A2.22,  may  be  rewritten  as 

T 

r  —  -t 

(A2.23)  C&(t)  =  'T  -  I  2  f(t)  f(t  +  t)  dt 

IT 

The  similarity  between  this  equation  and  the  computational 
formula,  equation  3.4*  is  evident. 


i  ->!•.  oc  /:  no  t.oco^b  jns*  sqq.6  ns  V  eb 
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,  N-L 

(3.4)  C(L)  =  ^  2i=i  G±  G1+l 

L  =  0,  1,  2,  (N-l)  . 

The  inability  to  use  infinitely  long  lags  means 
that  we  are  effectively  multiplying  the  true  autocorrelation 
(with  infinite  lag)  by  a  time  lag  window*  depicted  below. 


m 


m 


V 

:t) 

1 — 1 

II 

1  T' 

1  < 

T 

m 

V 

[*0 

1 

1  T  I 

1  = 

T 

m 

V 

(t) 

o 

H 

1  T  1 

1  > 

T 

m 

Taking  the  Fourier  transform  of  this  function  gives  the 
spectral  window. 


In  computing  the  power  spectrum*  the  result  is  thus  multiplied 
by  the  square  of  the  function*  Wq(^).  This  obviously  intro¬ 
duces  objectionable  side  bands  of  frequencies  into  the  compu¬ 
ted  power  spectrum*  as  the  first  side  band  has  an  amplitude  20 % 


of  that  of  the  main  lobe. 


a--jrc  .  :,'k  1  (J>° 


•If  (r,  , 


f  -)  0W 

ra 


.  V,  ■  L*.  1  J  l  '  I  3 
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To  partially  overcome  this  difficulty,  special  lag 

windows  which  give  less  weight  to  the  long  lags  have  been 

constructed.  These  are  treated  as  weighting  functions  for 

C  ,  producing  a  modified  apparent  autocorrelation  function, 
a 

(A2.24)  cma(T)  =  W(t)  Ca(T) 


Having  determined  this,  we  calculate  the  estimated  power 
spectrum  from  the  cosine  transform  A2.21. 


"[  QQ 

(A2.25)  P  =  —  f  C  (t)  cos  dx 

®  "  0  ma 

The  similarity  between  this  function  and  equation  3*5*  the 
computational  formula  is  also  evident: 


(3.5) 


N-l 

P(J)  =  2t  Z 


L=0 


W (l)  C(L) 


cos 


ttLJ 

M 


where  in  the  equation  t  is  the  digitizing  interval. 

The  estimated  power  spectrum  gives  a  weighted 
average  of  power  density  over  a  range  of  frequencies  and  not 
the  local  power  density.  As  the  maximum  lag,  f  ,  becomes 
greater  the  approximation  between  the  two  becomes  better. 


be  ''  c:  •  j  e  ;  o  ?of  '  ,  0 


( r£.£L) 


. 
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A2.3  The  Daniell  Window 

The  ideal  window  would  be  one  whose  frequency 
response  was  square-shaped: 


W(oi) 


W(o)) 

W(03) 


— ,  |  a)  |  <  L 

1  1  —  m 


2tt 


0,  |  “  |  >  - 
9  m 


The  Daniell  window,  used  in  the  power  spectral  calculations 
in  this  thesis,  approaches  this  ideal  as  closely  as  the  data 
allow.  To  calculate  its  time  lag  window,  we  calculate  the 
inverse  Fourier  transform  of  W(^),  above,  by  A2.l6. 


ir/m  •  /.w 

®(t)  =  ^  /  1  •  eltUT  doo 


-7r/m 


m 


27TiT 


r  i^T  -.Tr/m 
Le  ]  ' 

-Tr/m 


-  [e1  .  e-i">]  f  2i  (S) 


sin 


ttt 

‘m 


0)(T)  = 


7TT 

m 
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The  lag,  t,  can  have  values  between  0  and  -  T  ,  the  record 
length,  so  that  the  square  window  can  only  be  approximated. 

Use  of  this  Daniell  window  requires  more  computer  time  because 
the  autocorrelation  function  must  be  computed  for  all  lags 
up  to  the  maximum  record  length.  The  power  spectra  computed 
with  the  Daniell  window  represent  smoothed  averages  over  a 
controllable  range  of  frequencies  about  any  required  frequency. 


if  [■■■:  rfcr.i £  C.  wooed  ,i  - 
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APPENDIX  3.  LISTINGS  OF  COMPUTER  PROGRAMS 

This  appendix  documents  the  Fortran  IV  listings  of 
the  principal  computer  programs  which  have  been  developed  and 
used  in  this  thesis. 

A3.1  The  On-line  Printer  Plotting  Program 

This  program  gives  a  quick  method  of  plotting  (six 
traces  per  run)  on  the  on-line  IBM  1403  printer.  Input  data 
should  be  on  tape  with  the  reel  mounted  on  drive  number  4. 

The  control  cards  required  for  its  use  are  listed  below. 


1. 

123456X  -  beginning  in  column  1.  Columns  8  to  10 

inclusive  give  "RUNS",  the  number  of  graphs  to  be 

plotted,  including  those  that  are  skipped. 

2. 

Blank  card. 

3. 

Blank  card. 

4,  . 

5* . . •  1  card  for  each  run.  The  data  are  in 

(415)  format  and  gives  the  following  parameters: 

TOTAL  REJECT  CARDS  SKIP  where 

TOTAL  is  the  total  number  of  cards  in  the  run 

REJECT  is  the  number  of  cards  to  be  rejected 

before  the  plotting  begins 

CARDS  is  the  number  of  cards  to  be  plotted 

SKIP  is  the  frequency  of  plotted  points;  e.g.  a  ”2 

means  every  second  value  after  the  first  is  plotted 


sit  9riJ  8 Ui:‘  nl  beau 


-  :  .  .  :  i'  .  r.-en  j  x-v.)  Iof  .  )  ; r, 


v  s  BJJlonJt 

Ouori 

.  ::£  j  If  . 


T  I  .  JaTOT 


-  er{  3 n.i  j  !q  ari-t  eiclsa 
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$  JOB  7  3  5-21 

$  I  B JOB  GRAPH 

3IBFTC  PLOT  NODECK 

DIMENSION  IA ( 12- )  > IA2 ( 12 J)  >N ( 6  >20oJ )  > ISToRl ( 6 )  >  ISTOR2 ( 6 )  >MARK(6 ) 
DIMENSION  I  D  I  G I T ( 7 ) 

INTEGER  RUNS*TOTAL>REJECT >CARdS>CLASS>SKIm,S 
1  FORMAT ( 7A1 > I  3 ) 

3  FORMAT (415) 

4  FORMAT ( 8 - H 


15 

10 

98 

99 
102 


C 


r 

v. 


2 


C 

C 


C 

9 

C 


11 

C 

6 


201 

101 


1 


FORMAT ( 1 2  u  A 1 ) 


FORMAT ( 1X>18I4>7X) 


FORMAT  (  ]  2 OH]  *xxxxxx 

\L  SL  _y_  'L  X  -s  ->C  -><-  v_  M-  i  -H-  M-  J4 

|  /\  /\  A  A  A  A  A  TV  A  A  /V  A  A  A  j  7\  A  A 


X  X  J  -X- X  X  X  X  X  X  X  X  X  *  X  7\~  V\  -X"  "X'  -X-  X  7C 
a'  -X  -X  X  a-  X  XX  XX  X  -a-  X  X  X  X  J  X  X  X  X  A 


XXXXXXXXXXXX  X  X  ^  X X X X  X -X- X X  X  ) 


FORMAT ( IX  *  18 14 ) 


FORMAT (1UH1  CHANNEL  >14) 


FORMAT  ( 1 H I ) 

RE AD ( 5  > 1 )  (  IDIGIT ( 1 1 )  >  11  =  1 >7  )  >RUNS 
READS  CONFIGURATIONS  AND  NUMBER  OF  RUNS 
READ ( 5  >95 )  (  IA ( 17 )  >  17  =  1 > 12 J ) 

FOR  M AT ( 8  u  A  1 ) 

SETS  UP  BLANK  FIELu 


DO  2  11=1 » 12 u 
I A2  (  1 1  )  = I  A (  II) 

R  E  W I  N  D  4 

DO  5  11=1 > RUNS 

READ ( 5  >  3 ) TOTAL > REJECT , CARDS  > SKIP 
READ  ( 4  >  4 ) 

I  F  (  REJECT.  EG.TOTADGO  TO  55 
READS  PARAMETERS  FOR  THE  RUNS 

WRITE  ( 6  >  4 ) 

READS  HEADLR  CARD  FOR  L.ACH  RUN 
I R2= TOTAL-REJECT-CARDS 
C A RD S=CAR D S* 1 8 
DO  6  12=1 >6 

I F ( REJECT. EQ.O)GO  TO  9 
DO  7  I  3  =  1 >  REJECT 
READ  ( 4  >  8 ) 

REJECTS  THE  BEGINNING  OF  THE  CHANNEL 
RE AD ( 4  > 1 J )  ( N  (  I  2  *  I  3 )  >I3=1>CARDS) 

READS  SPECIFIED  NUMBER  OF  CARDS 
IF  (  I R2 • EQ • 0 ) GO  TO  6 
DO  11  I  3  = 1 > 1 R 2 
READ  ( 4  >  8 ) 

REJECTS  THE  END  OF  EACH  CHANNEL 
CONTINUE 


DO  101  I  T= 1 j  6 
WRITE(6>102) IT 

WRITE (6 >99)  ( N (  IT*  13)  >  13  =  1 > CARDS) 
DC  2  j 1  13=1, CARDS 
N (  I  T  >  I  3  )  =-N ( I T  > I  3 ) 

CONTINUE 


\ 


,  tr  )  T  A 

.  r  r  ,J.  III.  t  )  , 


(  :  )  TV  . 
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WR  I  TE ( 6  >98 ) 

C  LASS  =  5  0 
MARK ( 1 ) =1 1 
MARK ( 2 ) =31 
MARK ( 3 ) =51 
MARK ( 4 ) =7 1 
MARK ( 5 ) =91 
MARK ( 6 ) =1 1 1 

C  SETS  CENTERS  OF  GRAPH 

I NCR=CARDS~SKI P-1 
DO  12  I  4= 1 9  I  NCR  5  SK I P 
DO  13  15=1  ?  6 

13  ISTOR2 (  15 )  =  IDIGIT (  15 ) 

C  STOR2  HAS  THE  CONFIGURATIONS 

DO  14  16=196 
S  =  -  1 

IF (N (  16  9 14 ) .GE.G) S  =  1 

I F ( N (  16  9 14 )  . LT. 0 ) N ( 16  9 14 ) =-N (16  9  14) 

I STOR1 ( 16 ) =N (  I  6  9 14 ) /50 

I F (  IS  TORI (16). QT .9)1 STOR2 (I6)  =  IDIG1T(7) 
IF (  IS TORI (  16 ) .GE.9 ) ISTORl (  16 ) =9 
14  I S  TOR 1 (  I  6 ) =S  *ISTORl (  16 )+MARK(  16 ) 

Kl=  ISTORl (  1  ) 

K2=  ISTORl ( 2 ) 

K3=  ISTORl ( 3 ) 

K4=  I  ST  OR1 ( 4) 

K5=  ISTORl ( 5 ) 

K6=  ISTORl ( 6 ) 

I  A ( K1 )  =  ISTGR2 ( 1 ) 

I  A ( K2 )  =  ISTOR2 ( 2 ) 

I  A ( K3 )  =  ISTOR2 ( 3 ) 

I  A ( K4)  =  ISTOR2 ( 4) 

I  A ( K5 )  =  ISTOR2 ( 5 ) 

I  A ( K6 ) =IST0R2 ( 6 ) 
WRlTE(69l5)(lA(J5)9J5=l9l20) 

DO  16  18=1 9 120 
16  I  A (  I  8 )  =  I A  2 (  18) 

12  CONTINUE 
GO  TO  5 

55  CARDS=TOTAL*6 

DO  56  I  I  B B  = 1 > C A R D S 
56  READ  (498) 

5  CONTINUE 
STOP 
END 


$EI  TRY 
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A3. 2  The  Normalization  Program 

This  program  normalizes  all  traces  to  a  chosen  stan¬ 
dard  trace.  It  allows  for  input  traces  of  unequal  length  but 
the  output  traces  are  all  of  one  length.  A  formula  representing 
the  normalization  process  is  given  in  Chapter  3  as 

S.|G»td-0ral| 

(3'1)  "i  '  (°i  ’  51  Vo.  -S| 

The  control  cards  required  for  use  of  the  program  are  listed 
below  (symbols  are  defined  on  comment  cards). 

1.  IB,  STMC  in  (215)  format. 

2.  TESTHH,  STN08,  STN09,  CENTS T,  PTSST,  STCH,  SKI PI, 
NCARD,  MCI,  MC2,  MCN  in  (3(4X,Al),  1015)  format. 

3.  NRUN,  TESTH,  ZERO,  ONE,  TWO,  ...  ,  NINE  in 
(IX,  14,  11A1 )  format. 

4.  The  following  cards  are  required  for  each  run, 
giving  a  total  for  part  4  of  7  cards. 

a.  SKIP,  RECORD,  NOPTS,  MC,  MCC,  MN  in 
(715)  format. 

b.  CHAN,  CENPT  in  (215)  format  is 
required  for  each  of  6  channels. 


. 
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$  JOB  735o27 

$  IE  JOB  RMC8 
$  I BF TC  NORM  NODECK 

DIMENSION  VECT1 ( 5500 )  , VECT2 ( 5500 )  > NORM ( 6 )  ,T ITLl ( 5 )  >M( 6 ) 

REAL  VECT2  > NORM > MEAN  * M 

INTEGER  SKIP  >VECTl , RECORD , TOTAL 

INTEGER  R 3 33 3  >  STMC 

INTEGER  ADJUST  > CENPT > CHAN > SER I AL 

INTEGER  I  MAG (79) 

INTEGER  ONE>T/vO>THRLE>FOUR>FiVE>SiX>SEVLN>EIuHT>NlNE>ZLRo>TLSTFi 
REMIND  4 
REWIND  2 

DATA  C H 12/ 6  H  l  N  D  C  A  R  / 

READ  (  5  >  1 1 9  )  I  6  >  STMC 

C  IB .  J  IF  IT  IS  THE  FIRST  RUN 

C  1  IF  PART  Of-  OuIfUT  ThPl  HAS  BLcN  OBTAINED 

C  STMC  o  «.  STANDARD  NO*  OF  POINTS  10  THE  LEFT  OF  ClNTER  TO  bt  USlD 

119  FORMAT  (1515) 

101  FORMAT (A6) 

IF (  IB.EQ. 0 )GO  TO  100 
103  READ ( 2 > 101 ) CHARAC 

IF ( CHARAC. LQ.CH12 ) GO  TO  102 
GO  TO  103 

1 02  BACKSPACE  2 
100  CONTINUE 

INTEGER  CLuToT  >  A  d  J  >  P I S S I  >  S l CH>SKIP1 ?  o I NO  8  ?ST^^9  >TESTHH> T^T 
INTEGER  CHECjRECI 

DIMENSION  CHcC(8o)  >REC1 (550u)  >  REC2 ( 5  5  00 ) 

READ  ( 5  >110  )  TESTHH  >  ST N08  >  STN09  > CENT ST >  A  T  S  S I  >STCH>SKIP1 >  NC AKD  >  RiC  1 > 


11  0 


c 

c 

c 

c 

c 

Vw 

c 

c 

c 

c 

c 

c 


52 


7 


1MC2  >MCN 
FORMAT 
TESTHH 
STN08  . 
STN09  . 
C  l  N  T  SI 
FTSST  . 
S  T  C  H  .  » 
SKIP1  . 


(  3  (  4X > A 1  )  *  1 u I  5  ) 

.  .  AN  H 

«  .  THE  TENS  D I  'o  i  I"  OF  THl  STANDARD  RlOORD  nUMBlR 
**  THE  UNITS  DIGIT  OF  THl  SiANdARd  RlCORd  NUMbER 
*•  CENTER  OF  THE  STANDARD  CHANNEL 

•  «  NUMB t R  OF  Po I imT S  in  THE  STANDARD  CHANNEL 

•  •  NUMBER  OF  T Hi E  STANDARD  CHANNEL 

*•  Tile.  NUMBER  OF  RECORDS  THAT  MUST  Bl  f A Sold  dY 
THl. 


Tu 


VJ  L_ 


T  TO 


CORRECT 


STANDARD 


ReCoRD  N  U  iv1  B  l  R 


N C A R D  *  » •  THl  N U M B l R  OF 
MCI  *.  NO*  OF  POINTS  TO 
MC2  . .  NO.  OF  F0INTS  TO 
MCN  ,,  NO.  OF  POINTS  TO 
OBTAINING  THE 
DO  70  KZ=  1  >  99999  9 
READ  ( 4  ?  5 2  )  ( CHEC ( J )  > J= 1 > 8 0  ) 

FORMAT  ( 8  d  A1 ) 

IF  ( CHEC ( 2 ). lQ. TESTHH )  SKIP1=SKIP1 
IF  (SKIP1.LT. 0)  GOTO  71 
CONTINUE 


CARDS  IN  THl  STANDARD  RECORD 
lEFT  OF  Cl N TER  To  Bl  NORMALIZED 
RIGHT  OF  CENTER  TO  Bl  NORMAL  I Z lD 
EITHER  SIDE  OF’  C E N T E R  TO  lE  oSED 
NORMAL  1 Z 1  No  CONS T ANT 


IN 


71  IF  ( CHEC (8 ) .NE.STN08 )  WRITE  (6 >998)  CHEC ( 8 ) 
998  FORMAT  ( 1 X  > 14HCHEC ( 8 )  READS  >A1) 

IF  ( CHEC (9 ) .NE.STN09 )  WRITE  (6 >997)  CHEC (9) 
997  FORMAT  ( 1 X > 1 4HCHEC ( 9 )  READS  >A1) 

IF  (CHEC (8 ) .NE.STN08 )  GO  TO  69 


•  • .  I 


k  in  * >  l 


1  i  u  . . -  i 


( ...  .  ) 

3  A  (  IvX* 

H  . 

(  i  -  X  *  t  * )  - 

,  i 

•  •  1  •  T  uj 
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...  -t  i 

..S3 

'  «  5  )  -  i  (  .  '*) 

(  Li  l  T 

(  .  >  .  .  H  o  ) 

.1 

...  :  ) 


non 


1 3  2 


IF  ( CHEC ( 9 ) .  NE . STN09 )  Go  TO  69 
IF  (STCH.EQ.l)  GO  TO  59 
N3  =  NCARD* ( STCH-1 ) 

DO  5  8  KW  =  1 >  N  3 

58  READ  (4*53) 

53  FORMAT  ( 1 H 1 ) 

59  READ  (4*6d)  ( REC 1 ( J )  * J= 1  * P TSST ) 

6  -  FORMAT  (18  14) 

A  D  J  =  C  E  N  T  S  T  - 1  -M  C 1 
TOT  =MC1+MC2+1 
DO  61  KV= 1 ?  TOT 
IND=ADU+KV 

61  REC2 ( KV )=RLC1 ( IND ) 

A  D  D  =  0  • 

MCNN=2*MCN+1 
I S 1  =  MC1 +1 -MCN 
IS 2= IS1+MCNN-1 
DO  62  KU  = I  SI  *  I S2 

62  ADD  =  ADD+REC2 ( KU ) 

REAL  NORMAL  *  MEAN Y 
MEANY= ADD/FLOAT ( MCNN ) 

N  0  R  M  A  L  =  C  . 

DO  63  KT= I  SI > I S2 

63  NORMAL=NORMA L+ABS  (  Rc.C2  (  AT  )  -MEANY  ) 

REWIND  4 

RE  AD  (  5  *  1  )  NRUN  #  T  ESTH  »  Z ERu , ON  E  *  T  o  ?  TlIREl  ,  FuuR  ,  F  1  V  l  ?  S  I  X  *  SlVlN  ,  l_  i  onT  >  im 
1 1  NE 

1  FORMAT ( IX, 14* 11A1 ) 

NR  UN  ...»  THE  NUMBER  OF  RUNS 
T  E S  T H  • • .  AN  H 

ZERO... NINE  ARE  NUMBERS  o-9  INCL.  WITH  NO  SPACES 
SER I AL  =  0 
DO  2  KA=1 *  NR UN 
REWIND  3 

READ  (  5  *  3  )  SKI  P  ,  RlCORD  »  NOP  T  S  *  M C  >  MC C  *  PIN 

3  FORMAT (71 5) 

C  SKIP  ....  Nu.  OF  Ri_CORuS  TO  BE  SNiPPcD  AFTlR  PRlVIoUS  RECoRd  uur<L. 

C  RECORD  «.  THu  RlCORd  NUMBER  DcSIRld  i_X*  RECORD  3  IS  30 

C  NO RTS  ...  THu  NUMBER  OF  POINTS  IN  THIO  RtCGRD 

C  MC  . NO*  OF  POINTS  TO  THE  LEFT  OF  CENTER  T 0  dE  NORMALIZED 

C  MCC  *.  NO.  OF  POINTS  TO  RIGHT  OF  CENTER  TO  BE  NORMALIZED 

C  MN  ...  NO.  OF  PTS •  TO  EITHER  SIDE  OF  CENTlR  To  OBTAIN  CHAN.  NORM 

DO  4  KB  =  1  *  999999 
READ (4* 5) ( IMAG(o) *0=1,79) 

5  FORMAT ( IX  * 79A1 ) 

IF (  I M  A  G ( 1 )  • E  0 • 1 LSrH)SKIP=SKlP~l 
IF ( SKIP.LT . u ) GO  TO  7 

4  CONTINUE 

7  CALL  HE LP  (  Z L RO  , oim c: ,  T  Wo  ,  T i  iR E E  *  FoU R  *  F  i  V o  »  S  i  X  *  S u V l. * m  ,  E  I  ui  i T  ,  i  im l  ,  K l o o K u 

1  *  I  MAG) 

6  IF ( RECORD. GT. SERIAL) oO  TO  8 

14  wRITE(6>9) RECORD 

9  FORMAT ( IX, luHRECORD  NO. ,I4*15HIS  OUT  OF  ORDER) 

25  GO  TO  11 

8  SER I AL=REC0RD 
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DO  13  K E  = 1 *6 
RECORD* REC OR D+l 
READ (5 *3 ) CHAN,ClNPT 

CHAN  •  •  CHANNEL  NO*  — MUST  BE  WR 1  FlcN  WITH  RECORD  Nu*  In  FRONT 
CENPT  ...  THE  CENTER  POINT  OF  THE  CHANNEL 
I F ( RECORD • NE • CHAN ) GO  TO  14 
READ ( 4  5 16 )  (VECT1 (KM)  ,KM=1 *NOPTS) 

IF  (MC.EQ.STMC)  GO  TO  45 
LL=STMC-MC 
DO  46  J J= 1 ?  LL 

46  VECT2(JJ)=o. 

ADJUST =CENPT~1-MC 
TOTAL=MC+MCC+l+LL 
L J  =  L  L+ 1 

DO  47  JK=LJ, TOTAL 
JL= JK-LJ+1 
I NDEX=ADJUST+JL 

47  VECT2 ( JK)=VECT1 ( INDEX) 

GO  TO  48 

A  5  ADJUST  =  CENPT-1-MC 

16  FORMAT (1814) 

TOT  A  L=MC+MCC+ 1 
DO  17  KF= 1 >  TOTAL 
I NDEX=ADJUST+KF 

17  VECT2 ( KF) =VECT1 ( INDEX) 

C  VECT2  IS  FLGATINo  Pul  NT  VECTOR  VcCTl  IS  InTcGciR  VARIAdLc  VcCTuR 

43  S  U  M  =  C  • 

MNN=2*MN+1 
I S  T  AR  =  MC+ 1 -MN 
1ST AR2= ISTAR+MNN-1 
DO  18  KG= ISTAR * ISTAR2 

18  SUM=SUM+VECT2 ( KG) 

ME  A  N  =  S  UM  /  F  LOA  T  (  MN  N  ) 

S  U  M  =  0 

DO  19  KH= I  STAR  ? I  ST  AR2 

19  SUM=SUM+ABS ( VECT2 ( KH ) -MEAN ) 

M( KE ) =MEAN 

NORM ( KE ) =SUM 

WR I TE ( 3  > 1 2  )  (  I  MAG ( K I )  ,  K I  =  1  *  4  )  *  RECORD ,  (  I  MAG ( K I ) > K i  =  iO  *  7 9 ) 

12  FORMAT  ( IX  »  4A1 » I  5  *  7GA1  ) 

WR I TE( 3  *20  )  ( VECT2 ( KJ )  *KJ  =  1 *  TOTAL  ) 

2  0  FORMAT ( 1X*5E14. 7  ) 

13  CONTINUE 
REWIND  3 

DO  21  KK= 1 > 6 

READ (3* 31) ALPHA* RECORD? ( T 1 TLl ( KI )  , K I =1 , 5 ) 

31  FORMAT (A5> I5*5A6»4dH  ) 

READ (3 *20 ) ( VECT2 ( KJ )  *KJ  =  1 *  TOTAL) 

COR = NORMAL /NORM ( KK ) 

IF ( MC.EQ.STMC )  uO  TO  701 
LM  =  STMC-MC 
DO  702  JJ J  =  1 >  LM 
7^2  VECT 1 ( J JJ ) =  0 
LN  =  LM+1 
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DO  7^3  J J L  =  LN  »  TOTAL 

VECTl  (  JJL)  =  (  VECT2(  JJL)-M(KK)  )*COR 

7^3  IF  (  IABS  (  VECTl  (  JJL  )  ) . GT.999 )  VECT1(jjL)=999*(VECT1(JJL)/1ABS(VLCT11 
1  JJL )  )  ) 

GO  TO  24 

701  DO  22  KL= 1 >  TOTAL 

VlCTI  (  KL )  =  (  VECT2  (  KL)  -M  (  KK)  )  *C0R 

2  2  I  F  (  IABS  (VECTl  (  KL)  )  .  GT  .999  )  VECT  1  (  KL  )  =99  9*  (  VECTl  ( iCL  )/ IABS  (  VECT  1  (  nL  )  ) 
1  ) 

24  wR  I  TE  (  6  *  3  1  )  ALPHA  ,  RECORD  ,  (  T  I  T  l.  c  (  K I  )  ,  K  I  =  1 , 5  ) 

R3333=REC0RD-1 

I  F  (  K  K  •  E  Q  •  1  )w»RITE(2,31  )ALRHA,R3333  ,  (T  I  T  L  E  (  K  I  )  ,  K 1  =  1  ,  3  } 

WR  I  T  E ( 2  » 1 6  )  (  V  L  C  T 1  (  K  E  )  ,  KE=1, TOTAL) 

21  WR I TE( 6  »99 )  ( VECTl ( KE )  >KE=1 > TOTAL  ) 

99  FORMAT ( IX, 18 14) 

2  CONTINUE 

WRITE (2 ,101 )CH12 


GO  TO 

25 

1 1 

REW I ND 

4 

R  E  W  I  N  D 

o 

R  E  W I  N  D 

3 

69 

STOP 

END 

SIBFTC  SUB  DECK 

SUBROUTINE  HELP ( ZERO »u Nc» TWO, THReE, FOUR, P  IVc,SIX»SEVc.N,EIQHT  ,  N I  N  L  , 
1 R  E  C  0  R  D  ,  I M  A  G  ) 

INTEGER  AD, RECORD 

INTEGER  ONE,  r WO, THREE, TOUR, FIVE, SIX, SEVEN, ElGhT  ,  N 1 N  E , Z  E  R  0 
DIMENSION  1  MAG ( 79 )  ,N(5) 

RECORDED 
N( 1 ) =1^000 
N  (  2  )  =  1 J  0  J 
N  (  3  )  =  1 C  0 
N  (  4  )  =  1 0 
N  (  5  )  =  1 
DO  1  J= 1 ?  5 
K  =  4 + J 

I F ( IMAG(K) . EQ. ONE ) AD=1 
I F ( IMAG(K) .  EQ. TWO ) AD=2 
I F (  I  MAG ( K ) . EQ. THREE ) AD=3 
I  F  (  IMAG(K)  .  EQ.  FOUR  )*AD  =  4 
I  F  (  I  M  A  G  (  K  )  .  c_  Q  •  F  1' V  E  )  A  D  =  5 
I F (  I  MAG ( K ) . EQ.S1X)AD  =  6 
I  F (  I M A G ( K  )  • EQ • S l V E N ) A D= 7 
I F (  I  MAG ( K )  . EQ. E I GHT ) AD=8 
IT  (  I  M  A  u  (  K  )  .  E  Q  .  N  i  N  E  )  A  D  =  9 
I F ( IMAG(K) .EQ.ZcRO) AD=0 
1  RECORD=RECORD+AD*N  (  U  ) 

RETURN 

END 


SENTRY 


NORM 
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A3. 3  The  Autopower  Spectra  Program 

This  program  calculates  autocorrelation  coefficients 
and  autopower  density  spectra  of  input  data  with  variable 
format.  It  includes  a  choice  of  three  types  of  lag  windows: 
Hanning,  Hamming  or  Daniell.  The  relevant  formulae  are  given 
in  Chapter  3*  For  the  autocorrelation  coefficient. 


(3.4) 


C(L) 


1 


V'N  L  p  p 

Zi.  Ur.  Ur. 


N-L  1=1  i  i+L 


where 


L  =  0,  1,  2,  ...  ,  (N-l) .  The  formula  for  compu¬ 


tation  of  the  autopower  spectra  is 


(3.5) 


P(J)  =  2t  xjfl  W(L)  C(L)  cos  ^ 


where 


J  =  0,  1,  2,  .  .  .  ,  M 

M  is  the  specified  maximum  J  value 


Control  cards  required  for  the  use  of  this  program 
are  listed  below  (symbols  are  defined  on  comment  cards) . 

1.  INPUT,  NPOWER  in  (5X,  1015)  format. 

2.  (FMT(J),  J=l, 12)  in  (IX,  12A6)  format. 

3.  REC0R1,  REC0R2,  NOPTS,  MLAG,  Cl,  MTYPE ,  MR  in 
(5X,  1015)  format. 

DELT  in  (5X,  F10.5)  format. 


4. 


:  u.n  hsl'llot'%3  5  l;f  1 


. tf-tvic  (eioi  «} ;)  at  flawow  .twiti  ..x 

.4 


13  6 


$ JOB  735U59 

$  I B JOB  CROSS 


SIBFTC  ERK 1 

c  this  program  calculaTcS  autu  power  density  sPl_ctra  and  auto 
c  correlation  of  input  da  i  a  with  variable  formats 
COMMON/ DAT /A ( 3000 ) >B(3000)>C(6000)»D(6000) 

C  On  MON  /  C  Of  I  ST  /  R  ECOR 1  >  R  ECOR  2  >  NOP  TS  9  VILA  G  >  INPUT  ?  C  L  i  *  T  1  »  EHD 1  ?  C  2  >  C 1  >  iviT  YP  E 


DIMENSION  FMT ( 12  ) 

INTEGER  A  *  B 

INTEGER  RcCORl  >  R  E  C  OR  2  9  N  0  P  T  S  9  M  L  A  G  9  C  E  N  T  1  *  E  N  D 1  9  C  ^  »  C 1 
C  DELT  ***  THE  DIGITIZING  INTERVAL  IN  SECONDS  OR  MILES 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NOPTS  *** 
ML AG  *** 

END  1  *** 

CENTl  *** 
INPUT  *** 
FMT  DEFINES 
RECOR1  =  1? 
REC  OR 1  =  1 > 
C  2  *** 

*** 

Cl  ***  1 
*-**  3 

MTYPE  *** 
*** 
*** 

RATIO  NOPTS 
FOR  HANNING 


THE  NUMBER  OF  POINTS 
THE  MAXIMUM  FREQUENCY  INCREMENT 
THE  END  POSITION  OF  THE  CORRELATION 
THE  CENTRE  POSITION  uF  THE  CORK 


LAT I  ON 


VECTOR 
w  VECTOR 


UNIT  FROM  WHICH  DATA  IS  TO  BE  READ 
FORMAT  OF  INPUT  DATA  (1X,18Ia,7X) 
RECOR2  =  0  FuR  AUTo  CORRELATION 
R ECOR 2  =  2  FuR  CROSS  CORRELATION 
l  AUTO  CORRELATIONS  WANTED 
0  CROSS  CORRELATIONS  WANTED 
AUTO  POWER  DENSITY  SPECTRA  WANTED 
CROSS  POWER  DENSITY  SPECTRA  WANTED 

1  FOR  HA N N I N G  W I N D 0 W •  ( M L A G  SHOULD  E ( 

2  FOR  HAMM  I NG  WINDOW#  (MLAo  SHOULD 

3  FOR  DANIEL  W I ND  0 W »  ( M L A G  SHOULD 

/ M R  SPECIFIES  THE  IYPE  OF  D A N I c L 

AND  HAMMING  WINDOWS  MLAG  SHOULD 

TS 


i  \!  '■ 


N  P  0  W  E  R 


NO#  OF  DATA 


Ol 


PTS/10 ) 
EuUAL  NOPTS/ 10) 
EQUAL  i Ni 0 P T S  ) 

W  I  NDOW  (  RECOiviMEND  10  ) 
EQUAL  MR 


r  OR 


CALCUL, 


\ T I  0 1 m  OF  AOT OPuW'ER  S P E C T R 


t~\ 


1  FORMAT (  5X  > 10  I  5  ) 

8  FORMAT  ( 1 X  > 1 2 A6 ) 


9  FORMAT  ( 5 X,  FI  0.3) 

READ  (5*1)  INPUT  >  NPOWER 
READ ( 5  >  8 )  (FMT(U) >J  =  1 >12) 

READ  (5*1)  RECORl  >  R  E  C  0  R  2  >  NOP  TS  >  MLAG  >  C 1  >  i  -1 T  Y  P  L  >  M  R 
READ  (5 >9)  DELT 
CENT 1=MLAG+1 
END 1  =  MLAG+ CENT  1 
DO  5 1 * * * 5 * 7 J  I  K= 1 >NPOWER 
WRITE  (6*51)  IK 

51  FORMAT  ( 1H 1 > 4HSET  ?I2>14H  OF  INPUT  DATA) 

R  E  A  D (  INPUT j  F  M  T )  ( A ( U )  >U=1> NOPTS) 

WRITE  ( 6  >  F MT )  ( A ( U )  >  J  =  1 ?  NOP T S ) 

C  2  =  0 


I F ( R tCOR2 • NE • 0 ) GO  TO  5 
DO  6  0=1?  NOP  T S 
6  B  (  U  )  =  A  (  U  ) 

C2  =  1 
GO  TO  7 


5  READ ( INPUT  > FMT )  ( B ( U )  > J=1 > NOPTS ) 

WRITE  (6* FMT)  (3(0) >J=1> NOPTS) 

7  CALL  CROSS ( CENTl » NOPTS  >  C2  > MLAG  »M TYPE ) 

I  F  (  C2  .  tQ.  0  )  CALL  OUTPUT  (  M lAo  >  E N D 1  >  RECORl  9  R l. C u R 2  9  0  ) 

IF  (  C2  •  EQ  •  1  )  C  A  L  L  OUTPUT  ( i’,i  l.  A  u  >  _  1'fDl  9  \ \ E  C  Q  R 1  9Rt_CuR2  9iviLAu) 
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DO  2 0  J  A  =  1  * M L A G 
JB  =  CENT 1+JA-1 

20  WRITE  (3*21)  JA»C(JB) 

21  FORMAT  (IX, 15, 2X, El 4. 7) 

WRITE  (3*22) 

22  FORMAT  (6H  END  ) 

FN=  1 • / ( 2 • *DL  LT ) 

WRITE  (6*10)  DELT  *  FN 

10  FORMAT ( IX , 5HDELT= ,F8 . 3  *  1 X , 1 8 HNYQU I  ST  FREQUENCY =  ,F8 . 3 ) 
WRITE  (6*12)  NOPT  S  >  MLAG  »MR 

12  FORMAT ( IX  , 6HNOPT  S= » I  5 , 1 X  ,  5  H  M  L  A  G  =  , 15, 1X,3HMR=, 15 ) 

CALL  POWER (MLAG, DELT ,CENT1 , M  TYPE  » MR ) 

DO  11  J  =  1 ,  END1 

11  C  (  J  )  =  D  (  J  ) 

I F ( Cl . EQ. 1 ) CALL  uUTPUT ( MLAo  ,END1 ,MTYPE ,C1 , MLAG ) 

DO  23  JA= 1 > MLAG 
JB  =  CENT  1+ JA-1 

23  WRITE  (3,21)  JA,C(J6) 

WRITE  (3,22) 

50  CONTINUE 
STOP 
END 


$  I 3FTC  SU31 j 

SUBROUT INE  CROSS(CENTl , NOP TS  >  C2  > MLAG ,MTYPE ) 

COMMON/ DAT /A (3000) ,B(3000),C(6000),D(600Q) 

INTEGER  A , B 

INTEGER  C E N T 1 , N 0 P T S , C 2  > M L A G 
KX  =  0 
KY  =  0 

DO  3  J= 1 , NOPTS 
KX  =  KX  +  A ( J ) 

3  KY  =  KY  +  B ( J ) 

KX  =  KX/NOPTS 
KY  =  KY/NOPTS 
DO  4  J=l, NOPTS 

A ( J  )  =  A ( J )  -  KX 

4  B ( J  )  =  B( J )  -  KY 
WRITE  (6,6) 

6  FORMAT  (1HL»38HINPUT  VALUES  WITH  AVcRAul  v • C  *  R  E M 0  V  E  u ) 
WRITE  (6,5)  ( A ( J ) *  J=l, NOPTS) 

5  FORMAT  (IX, 1814, 7X) 

CALL  DOT ( NOPTS , 0  ,  C ( CENT  1 )  , 1 , M T Y P E ) 

DO  1  J= 1 »  MLAG 
K1 =CENT 1+ J 

CALL  DOT ( NOPTS , J ,C ( Kl ) , 2  * M T  Y P E ) 

I F ( C2.NE. 0 ) GO  TO  1 
K2=CENT1-J 

CALL  DOT (NOPTS, J ,C ( K2 ) ,3 ,MTYPE ) 

1  CONTINUE 
RETURN 
END 
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$  I  B  c  T  C  SUB  1 1 

SUB  ROUT  I  j  -4  E  D  0  T  (N*NS>ANS> MODE? MTYPE) 

COMMON/ DAT /A ( 30 CO ) ,6(3000) >C(6000)»D(6000) 

INTEGER  A, 6 
F  N  U  M  =  0  . 

L  =  N  -  N  S 
DO  1  U= 1 »  L 
K  =  U  +  N  S 

1  FNUM  =  FNUM  +  FLOAT (A(U) )* FLOAT (B(K) ) 

I F ( MTYPE. EQ. 3 ) GO  TO  2 
AN  S  =  FNUM/ ( FLOAT ( L )  ) 

GO  TO  3 

2  AIMS  =  FNUM/ (FLOAT  (N)  ) 

3  RETURN 
END 

SI3FTC  SU814 

SUBROUTINE  POWER  (  MLAG  , DELT  >CEN'T1  ,MTYPE*MR) 

COMMON/ DAT /A ( 3C0G ) ,B(3000) , C(6000) ,D(6000) 

INTEGER  A  >  B 
INTEGER  CENT  1 ,END1 
DO  1  L= 1 » MLAG 
L 1  =  CENT  1  +  L 

CALL  W I N  D  0 w ( M R  ?  L  > W ?  M I Y  P  E ) 

1  C(L1)=C(L1)*W 

END  1  =  MLAG  +  Cl.  NT  1 

CALL  OUTPUT (MLAG *  END 1 ,1 >  0  >  ML AG ) 

DC  =  G  . 

DO  2  L= 1 » MLAG 
L 1  =  CENT 1  +  L 

2  DC  =  UC  +  C  (  Ll  ) 

D(CENTl)  =  DELT*C ( CENT1 )  +  2.*DELT*DC 
DO  3  J  = 1 > M L A G 
A  C  =  0  • 

DO  4  L  =  1 » MLAG 
Ll  =  CENT  1  +  L 

CARG1  =  3 .  141 59* ( FLOAT ( J )  ) * ( FLOAT ( L )  ) / ( FLOAT ( MLAG )  ) 
4  AC  =  AC  +  C ( Ll ) *COS ( CARG1 ) 

J1  =  CENT  1  +  U 

3  D(U1)  =  DELT*C(CENTl)  +  2**DtLT*AC 
RETURN 

ENp 


$  I  E  FTC  SUB15 

SUBROUTINE  WINDOW (MR * L,W> MTYPE ) 

CARG  =  3. 14159* ( FLOAT ( L ))/( FLOAT (MR ) ) 
GO  TO  (5*697) » MTYPE 

5  W  =  0.5  +  0. 5*COS ( CARG) 

GO  TO  8 

6  W  =  _.54  +0.46*COS(CARG) 

GO  TO  8 

7  w  =  ( SIN ( CARG) ) /CARG 

8  RETURN 
END 
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T  '  (  £  .  .  iY  T  )  i  I 

(  ( .  )  I 

_j  *  ■  i  i 

)  it  {  )D.(  )  .  (  ..  )  \T- 

2«A 

_  *r*_ 

■ .  T  t  <  if 

•  . .  =  _  J 

.  =OA 

{  ’  V  / 

I  I  f  J  f  ) 

)  ..  H)  \  (.(  )  T  : )  X .  ■  *  0 

r  e  >  ♦  ) 

)  -r  -  V  •  +  •  = 


139 


SIBFTC  SUB  1 2 

o^BRuUT  1 1 'l £  OUTPUT  (iv'ILAG»c:i'iDl  >  R E C U R 1  *  R  E  C  U  R  2  »  i  A u  ) 

COMMON/ DAT /A ( 3000 ) ,  B(30  00) ?C(6000) *D(6Q00) 

INTEGER  SKIP >END 1 , RECUR 1  *  RECUR 2 

4  FORMAT ( 6 ( IX* l6?lX>lH*E10»3>lX»2ri  )) 

5  FORMAT (6(8X>1H*11X»2H  )) 

CALL  PAGE ( L7  *  SKIP , RECUR  1 *RECGR2 ) 

I R=-MLAG+ I  AD 
Kl=IR-6 
K  2  =  I  R  -  5 
K3  =  I  R-4 
K4  = I R- 3 
K3  =  I R-2 
K6=  I  R-l 
J= 1  +  IAD 
12  Kl=Kl+6 
K2  =  K2  +  6 
K3=  K3  +  6 
K4=  K4+6 
K5  =  K5  +  6 
K6  =  K6  +  6 
L  =  J 

Ll  = J+l 
L  2  =  J  +  2 
L  3  =  J  +  3 
L4  =  J  +  4 
L5  =  J  +  5 
J=  J  +  6 

I F ( J.GT.END1 )  GO  TO  11 

WRlTE(6>4)Kl  >  C  (  L  )  »  K  2  ?  4  (  E 1  )  ?  K  3  *  C  (  L  2  )  »  K  4  *  C  (  2  3  )  ?  4  5  *  C  (  L  4  )  *  K  6  *  v_  (  L  5  ) 

S K  I  P  =  SK I P+ 1 

L7=L7+1 

IF  (L7.EQ.45)  CALL  PAGE ( L 7 , SK I P , RECUR  1 , RLCOR2 ) 

IF  (SKIP.EQ.5)  WRITE (6*5) 

I  IF  (SKIP.EQ.5)  SK I P  =  4 
GO  TO  12 

II  I R  =  MOD ( E N D 1  + 1 - L  >  6 ) 

IF  (  I  R  •  EQ  .0)  w  R I  T  L  (  6  *  4  )  4 1  *  C  (  E  )?  K2  >  C  (  L 1  )  *  K 3  >  C  (  L )  *  K4  ?  C  (  L  3  )>  4  5  * 
1  C ( L4) *K6>C( L5 ) 

IF  (IR.EQ.O)  RETURN 
GO  TO  (6*7*8>99lG) » I R 

1  w  W  R  I  T  E  (  6  » 4  )  Kl *  C ( L )  *  K2  *C ( Ll  )  »  43  >C ( )  ,  \4  *  C  (  l-3  )  »  45  >  C ( L4 ) 

RETURN 

9  WRITE(6*4)K1*C(L) ,K2»C(L1) , K3 * C ( L2 ) , K4 , C ( L3 ) 

RETURN 

8  WRITE(6»4)K1*C(L)»K2*C(L1)*K3*C(L2) 

RETURN 

7  WR I TE ( 6  *4) Kl >C ( L )  *K2  »C ( Ll  ) 

RETURN 

6  WR  I  TE  (  6  *4  )  K1  *-C  (  L  ) 

RETURN 

END 
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SIBFTC  SU613 

SUBROUTINE  RAGE  (LI  ,  Six  I  R  ,  RLCuRl  ,RECuR2  ) 

INTEGER  RECORl * RECOR2 

ISKIP=G 


1 


4 


7 

8 
1 

2 


L 1  =0 

IF  ( R  E  C  0  R  2 • E  G • j )  W  R  I  T  E  (  6  » 1  ) RECUR 1 ' ,  R  E  C  0  R 1 
IF  ( RECOR2 - EQ. 2 )  WR I T  E ( 6  *  1 ) RECORl ♦RECURS 

FORMAT  (1H1,29X>28HC  ROSS  CORRELATIONS  uF  RECORd,1X»I5,1X»11H  V;  I  T  h  R  t 
1  CORD » IX* I  5 ) 

IF  ( RECOR2.EQ. 1 )  WRITE (6 >4) 

IF  ( RECOR2.EQ. 3 )  WRITE! 6 >5) 

FORMAT  (1H1>29X»26HAUT0  R  u  W  E  R  D  E  N  S  1 1  Y  SRlCIRA) 

FORMAT  (  1  FT  1  >  29X  »  27HCR0SS  ROWER  DENSITY  SPECTRA) 


IF(RECOR2»i_QeO)uO  TO  6 
IF(RECOR2.EQ*2) GO  TO  6 
IF  ( R E COR  1 • EQ • 1 ) WR I T  E ( 6  » 7 ) 

IF  ( RECORl. EQ. 2 )  WRITE  (6,8) 

I  F  (  RECOR  1  o  EQ .  3  )  w'R  I  T  E  (6,9) 

FORMAT ( 30X *19HW I TH  HANNING  WINDOW) 

FORMAT  (  30X  *  1  9HV/I  TH  HAMMING  WUNDOW  ) 

FORMAT ( 3 Q X  ? 1 8 H W I T H  DANIEL  WINDOW) 

WRITE(6>2) 

FORMAT  (  3  Q  X  5  17  3  H  ***********  ^  ^  *  *  *  *  ''r  ^  ''  ^  Sc  7*r  tc  Tf  7C  75"  i<r  75"  7 r 

1  ///) 


~/<r  Sc  Sc  Sr  Sc  Sc  “X“  "X*  Sc  7V1  Sc  Sc  Sc  Sc  7 S 


WR I TE ( 6  » 1 0 ) 

1C  FORMAT (6(6X,1 HU, 8X,4HP( U ) ,3X) ) 
GO  TO  11 
6  WR  I  T  E  (  6  >  2  ) 

W  R  I  T  E  (  6  >  3  ) 

3  FORMAT ( 6 ( 6X»1HL,8X,4HC( L)  , 3  X )  ) 
11  RETURN 


END 


SEN  TRY 
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A3. 4  The  General  Stacking  Program 


This  program  enables  the  compositing  of  a  maximum 
of  8  traces  per  stack  with  variable  weighting  factors.  The 
sub-routine  "SUB2",  called  "CODER" *  is  useful  as  it  enables 
two  numbers  in  14  format  to  be  stored  as  one  number,  thereby 
doubling  core  storage  for  the  input  data.  Control  cards  required 
for  its  use  are  listed  below  (symbols  are  defined  on  comment 
cards ) . 

1.  NCAKDS ,  NPTS,  NRUNS,  NSKIP,  INDEX,  INPUT,  IPTST, 
IPTEND,  JOUT  in  (5X,  1515)  format. 

2.  The  next  control  cards  depend  on  the  parameter 
INDEX. 

a.  INDEX  =  1.  For  each  run  three  cards 
are  required. 

i.  NST  in  (5X,  1515)  format, 

ii.  IRWD,  (NTRACE(I),  1=1, NST)  in 
(5X,  1515)  format. 

iii .  (WT(I),  1=1,  NST)  in  (5X,10F6.5) 

format . 

b.  INDEX  =  2.  Two  cards  are  required. 

i.  NST,  NS TART  in  (5X,  1515)  format, 
ii.  (WT(I),I=1,NST)  in  (5X,10F6.5) 


format . 


:a  .i«=i  ,  [ 
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98 

97 

96 


95 


C 

C 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IT  IT 


o 


$  JOB  7  35^43 

$  IE  JOB  R M C 9 
$  I B FTC  GENSTK  NODECK 

99  FORMAT  (  1  H  1  *  2  3  H  THIS  IS  STACK  NO.  »I2»50H. 

1PERCENTAGES  OF  THE  VALUcS  OF  ?12*10H  CHANNELS • ) 

FORMAT  ( 1H1,4X,4HX( I ) ,5X,4HStT  *1L*5X*4HSlT  *  i  2 , 5  X  ,  4piSlT 
FORMAT  (1HK) 

FORMAT  ( 8 0H1 . . . SIX  SlTS  dP  STaCnlu  DATA  FuLLuw 

1 T ATEMENT •  . .  .  ) 

FORMAT  (8  OH  000 0 0000000 00 0  000 00 0000 00000  0 000000 00 00000000000000000 
1000000000000000000  ) 


.OimPuSED  oF 


) 


T  h  i  S 


1 

2 

3 

4 

5 

6 

7 

8 
9 


1 


FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 

FORMAT 


(  5X  *  1 5  I  5 ) 

( 5X* 10F6. 5 ) 

(  8 OH 

( IX* 18I4*7X ) 

(IX* I4*1X*F4.0) 

( 6H  END  ) 

(5X* 14,3 ( 5X* 16  )  ) 


(4(  IX* 14  )  ) 

(  2  (  3  X  >  I  7  )  ) 

DIMENSION  W  T (10)  ,  TLMPOR ( 2 *10)  *  NT RACE ( 10 ) 

COMMON  A  (  8  *  2  0  0  0  ) 

INTEGER  A*  TEMP ( 2  *  1 0  )  , B ( 2  *  2 000 )  *  GRAPH ( 1 0 )  *  RESULT  ( 400 0 ) * HuLu ( 20 ) 
EQUIVALENCE  ( NT RACE ( 10)  *  GRAPH ( 10)  )  *  (3(2*2000)  > RESULT ( 4000 )  )  > ( TEMP ( 
12*10)  *  HOLD ( 2 J )  ) 


READ  (5*1) 
NC ARDS 
NPTS 
NRUNS 
NS  T 
INDEX 
I  N  D  E  X 
N  S  K I  P 


NCARDS, NPTS, NRUNS  >NSKIP  ,  INDEX*  INPUT  *  it^TST  *  IPTEND»^OUT 


1 

1 


=  ^ 


NUMBER 
NUMBER 
N  U  M  i3  E  R 
NUMBER 
STASK  INST 
STACK  IN  S  T 
N  U  M  8  E  R 


OF  CARDS  IN  ONE  TRACl 
OF  POINTS  IN  ONE  TRACL 
OF  RUNS  TO  uE  MADE 
OF  TRACtS  TO  BE  USuD  PER  STaCK 
TRACES  TO  OBTAIN  dNE  TRACl 
TRACES 


IPTST 
I PTEND 
JOUT 
JOUT 
NREAD  = 
I  F 

N  WR  I  T 


u  F  TRACES 
©FIRST- NST/2@ 
POINT  NUMBER  TO 
POINT  NUMBER  AT 

1  USE  AUTOPLOT 

2  USE  PLOTTING 


SIVElY  Tu  OBTAIN  (nST-1)  TRAC; 
TO  dc  SKiPPlU  blFuRl  AKRiVinG  aT 
TRACE 

AT  FOR  pLuT 
TO  End  PLOT 


.D 


START 


OUTPUT 


( MOD (NREAD, 6 ) . EQ. 


0 


WHICH 
AS  OU 

PROGRAM  735021  AS  OUTPUT 
0)  READ  (INPUT, 3) 


WRITE  (2*96) 

IF  (NSKIP.EQ.O)  GO  TO  12 
DO  11  I A  = 1 *  N  S  K  I  P 
DO  10  I  B  = 1  * NCARDS 

10  READ  (INPUT, 3) 

NREAD  =  NREAD  +  1 

IF  ( MOD ( N  R  E A  D , 6 ) • EQ • 0 )  READ  (INPUT, 3) 

11  CONTINUE 

12  GO  TO  ( 30  *  5w )  , INDEX 
30  DO  29  NR=1, NRUNS 

READ  (5*1)  NST 


READ  (5*1)  IRWD* (NTRAC£( I ) » 1  =  1  ,  NST ) 
READ  (5*2)  ( w  T ( I )  , I  =  1 , N  S  T  ) 

IF  (  I R ft D • E Q • 1 0 )  RE l I N D  4 
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IF  (IRWD.EQ.10)  NREA D  =  C 

C  IRwD  =  1j  TAPE  MUST  3E  REWOUND  BEFORE  ANY  i'-ioK£  STACKING  C  A  in  OCCUR 

NPTS2  =  NPTS/2 
DO  16  1=1 >  NST 
K A  =  N  TRACE (  I ) -NREAD-1 
IF  (KA.EQ.o)  GO  TO  13 
DO  15  I A= 1 ?  K  A 
DO  14  I B  =  1  » N  CARDS 
14  READ  (INPUT, 3) 

NR  E AD  =  NREAD  +  1 

IF  ( MOD (NREAD, 6 ) • EG. J  )  READ  (INHUT, 3) 

35  CONTINUE 

13  CALL  FINDER  ( I , NCARD5 ,NPTS , INPUT ) 

NREAD  =  NREAD  +  1 

IF  ( MOD (NREAD, 6 ) • LQ. 0  )  READ  (INPUT, 3) 

16  CONTINUE 

C  T H E  @NST@  SETS  OF  DATA  ARc  NOW  STORED  IN  i v i A T R  1  X  A 

DO  18  J=1 , NPTS2 
JA=2*U-1 
U  B  =  2  *  U 
RES  1  =  C. 

RES  2  =  0. 

DO  17  1=1, NST 

CALL  CODER  ( A (  I  ,  J )  , TEMP (1,1)  ,  TEMP (2,1)  ,  1  ) 

TEMPOR (1,1)  =  WT ( I )*FLOAT( TEMPI  1 , I )  ) 

TEMPOR ( 2 » I )  =  WT ( I ) *FLOAT ( TEMP (2 , 1 )  ) 

RES1  =  RES  1  +  TEMPOR (1,1) 

17  RES2  =  RES 2  +  TEMPOR (2, I) 

RESULT (JA)  =  I  NT (REST ) 

RESULT ( UB )  =  I  NT ( R ES2 ) 

18  CONTINUE 

C  THE  VECTOR  gRESUl-TG  NO1/;  CONTAINS  THE  STACKING  OF  i n S T  CHANNELS 

NF  =  2*NP  T  S2 
WRITE  (6,99)  NR, NST 
WRITE  (6,4)  (RESULT (  I  ),  1  =  1  ,NF  ) 

GO  TO  ( 24,25 ) ,UOUT 

24  WRITE  ( 3  >  9 )  ( I , RESULT ( I ) , I  =  1 , NF ) 

WRITE  (3,6) 

GO  TO  29 

25  WRITE  (2,4)  ( RESULT (  I  )  , I  =  1 , NF ) 

N  W  R  I  T  =  N  W  R  I  T  +  1 

IF  ( MOD ( N w R I T  » 6 ) • E Q • 0  )  WRITE  (2,96) 

29  CONTINUE 

N ZERO  =  6  -  MOD ( NR UN S, 6) 

IF  (NZERO.EQ.  6)  GO  TO  28 
DO  27  1=1, NZERO 
DO  26  11=1* NCARDS 

26  WRITE  (2*95) 

27  CONTINUE 

28  GO  TO  80 

55  READ  (5,1)  NST , NST ART 

C  NS  TART  THE  NUMBER  Tu  START  STACKING  FROM  AFTER  SKIhhINu 

C  @NSK I P@  RECORDS 

READ ( 5  > 2 )  ( WT ( I )  ,  1  =  1  ,NST ) 

NPTS2= NPTS/2 
N  S  T  2  =  N  S  T  /  2  +  1 
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DO  49  NR  =  1  > NRUNS 

IF  ( (NR.EQ.l ) .AND. ( NS T2 . EQ . NST AR T ) )  GO  TO  34 
IF  (NR.GT.l)  GO  TO  36 
DO  31  1=1 9  NST 

IF  (  ( I .GT. 1 ) .AND.  (  1  • LE.NST2 )  )  GO  T o  32 
CALL  FINDER  (  I  9  INC  A  R  D  S  >  N  F  T  S  ?  INPUT  ) 

NREAD  =  NRcAD  +  1 

IF  ( MOD ( NREAD 9 6 ) • EG  •  0 )  RcAD  (INPUT >3) 

GO  TO  31 

3  2  DO  3  3  I A= 1 ’NPTS2 
33  A (  I  *  I A )  =  A ( 1 ? I  A ) 

31  CONTINUE 
GO  TO  40 

l  '+  DO  35  1=1  9  NST 

CALL  FINDER  (  I  9NCARDS,iMFTS9  IimPUT  ) 

NREAD  =  NREAD  +  1 

IF  ( MOD ( N  R  E  A  D  9  6 ) • E  Q • 0 )  READ  ( INPUT  93 ) 

35  CONTINUE 
GO  TO  40 

3  6  NST  L 1  =  NST-1 
DO  38  1=1 9NSTLI 
DO  37  I  A  = 1 9  NP TS 2 
I  1=1  +  1 

37  A (  I  9  I  A  )  =  A( I  I  9  I  A ) 

38  CONTINUE 

IF  ( (NRUNS-NR) .LT, (NST2-1 ) )  Go  To  39 
CALL  FINDER  ( NST 9  NC A RDS  9 NP T S  9  INPUT ) 

NREAD  =  NREAD  +  1 

IF  (  MOD  (  NRt_AD  9  6  )  •  EG  •  g  )  RcAD  l  INPUT  93) 

GO  TO  40 

39  DO  43  I  A  = 1 9NPTS2 

43  A ( NST  9  I  A )  =  A { NSTLl 9  I  A ) 

C  THE  DATA  IS  NOW  IN  CORE  FOR  1  RUN 

40  DO  41  J  =  1 » NPTS2 
JB=2*J-1 

JC  =  2* J 
POINTl=o. 

POI NT2=0. 

DO  42  J  A  = 1 9  N  S  T 

CALL  CODER  (A(JA9J)  9  1 E  M  P  (  1  9  jA  )  9  IEIviP(29JA)  9I  ) 

T  E  MPOR  (  1  9  j  A  )  =  FLOAT  (  TEMPI  1  9  j  A  )  )  »•«  J  (  j  A  ) 

T  E  M  P  0  R  (  2  9  J  A  )  =  FLOAT  (  F  E  PiP  (  2  9  J  A  )  )  *  /*  T  (  J  A  ) 

POINT1  =  Pol  NT 1  +  TEMPOR(l9JA) 

42  POINT2  =  PO I  NT  2  +  TEMPOR(2  9JA) 

IPT1  =  INT(POINTI) 

I  P  T  2  =  I N  T  (  P  0  I  N  T  2  ) 

IM  =  MOD ( NR  9  3 )  +  1 
GO  TO  ( 44945946 ) 9 IM 

44  CALL  CODER  ( A ( 8 9 J )  9 1 PT  1  9  I pT 2 9 2 ) 

GO  TO  41 

45  CALL  CODER  ( B ( 1 9 J )  9 1 PT 1 9  I pT 2 9 2 ) 

GO  TO  41 

^3  CALL  CODER  (  uj  (  2  9  J  )  9  1  P  T  1  9  I  P  T  2  9  2  ) 

41  CONTINUE 

AFTER  ANY  3  RUNS  THERE  wILL  BE  0  OEFS  oF  STaCnlD  DATA  STOKtiD  Im 

A ( 8  9  ) 96(19 )  98(29  )• 
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I F  ( MOD ( NR  *  3 ) . NL . 0 )  uO  TO  49 
60  TO  (51*52) * JOUT 

51  N  2  =  N  R  -  2 
N 1  =  N  R  - 1 

WRITE  (6*98)  N2  *  N 1 9  NR 
WRITE  (6*97) 

DO  47  J=1 9  NP  TS2 
JB=2*J-1 
J  C  =  2  *  J 

CALL  CODER  ( B ( 1 > J ), GRAPH ( 1 )» GRAPH ( 4 )» 1 ) 

CALL  CODER  ( B ( 2 » J )  >GRAPH ( 2 )  *  GRAPH ( 5 )  > 1 ) 

CALL  CODER  ( A ( 8 * J )  > GRAPH ( 3 )  *  GRAPH ( 6 )  * 1 ) 

WRITE  (6*7)  JB  *  GRAPH ( 1 )  » G  R  A  P  H  (  2  )  >  G  R  A  P  H  (  3  ) 
WRITE  (3*8)  JB  *  GRA PH ( 1 )  *  J RA  P H ( 2 )  *  G RA  PH ( 3 ) 
WRITE  (6*7)  JC  >  GRAPH ( 4 )  > GRAPH ( 5 )  » GRAPH ( 6 ) 
WRITE  (3*8)  JC* GRAPH (4) > GRAPH ( 3 ) >gRAPH(6 ) 

47  CONTINUE 

WRITE  (3*6) 

GO  TO  49 

52  DO  55  K  =  1  *  2 

DO  54  J  =  1  * N P T S 2 

CALL  CODER  ( B ( K * J )  * G RA PH ( 1  )  > GRAPH ( 2 )  *  1 ) 

MAK  =  MOD ( J  *  9 ) 

IF  (MAK.EQ.O)  MAK =9 
JB  =  2* MAK  -  1 
JC  =  2*MAK 
HOLD(Jb)  =  GRAPH ( 1 ) 

HOLD ( JC )  =  GRAPH ( 2 ) 

IF  ( ( (NPTS2-J) . LT.9 ) .AND. ( JC.NE. 18 ) )  oG  TO  59 
IF  (JC.NE. 18)  GO  TO  54 
GO  TO  61 

59  IF  (J.NE.NPTS2)  GO  TO  54 

6 1  WR  I  T E  (6*4)  ( HOLD ( JM )  > JM= 1 , JC ) 

WRITE  (2*4)  (HOLD(JM)  *  J  Fi  =1  *  J  C  ) 

54  CONTINUE 

NWRIT  =  N  WRIT  +  1 

I F  ( MOD ( N  W  R I  I  ?  6 )  . E G . 0 )  W RITE  (2*96) 

WRITE  (6*97) 

55  CONTINUE 

DO  56  J  =  1 *  NPTS2 

CALL  CODER  ( A ( 8 > J ) >GRAPH ( 1 ) , GRAPH ( 2 ) » 1 ) 

MAK  =  MOD ( J  9  9 ) 

IF  (MAK.EG.U)  MAK=  9 

JB  =  2*MA K  -  1 

JC  =  2*MAK 

HOLD ( JB )  =  GRAPH ( 1 ) 

HOLD ( JC )  =  GRAPH ( 2 ) 

IF  ( ( (NPTS2-J) .LT.9 ) .AND. ( JC.NL. 18 ) )  ou  TO  60 
IF  (JC.Nu.18)  GO  TO  56 
GO  TO  62 

6 d  IF  ( J . N E • N P T S 2 )  GO  TO  56 

62  WRITE  (6*4)  (HOLD(JM)  ,Jivi=l,JC) 

WRITE  (2*4)  ( HOLD ( JM )  *  JM= 1 , JC ) 

56  CONTINUE 

NWRIT  =  NWRIT  +  1 

IF  (MOD  (NWRIT  96)  .EG).  0  )  WRITE  (2*96) 


(  •  •  t  *  .  )  J  ) 

I  »  u  ^ 

•  ,U«  •  )  .-  U  «  (  w  «  )  ) 

«  (  S  )  «  v  L  f  )  )  !  j  I 

-  ♦  (  )  'V  »  (  )  / 

-y  ♦  l  i  1  *“»  .  t  L  v  ’  ♦  l  )  J  >  ?1  . 

♦  t  .  )  ^ 

„  f  l  r1  )  :  .  “< 

( d  *  )  ;  i  :" 

S  «I*> 

f.  [  ’  i  : ,  i  c  ■«  }  i 

(e.LfUOi  =  y  A 

,  1.  (  .  -,.JL) 

£i  (s-th,  .:u;.u 

{  DLf  ;  *r,L«  ( 

-H  .  )  (♦>«£) 

1  •  k  •  l  5  •'  i 

(  .  e*  3)  I 

)  « ( )a)  . .  j  J-  O 

(  ♦  L  ) 

(  .  .  .  ) 

(  )  4  /•  /.  )  =  (  L )  J 
(  , )  r  /,>!£  =  (  3^  )  0  J 

i(  •  <  .  J  »  .  •.  •  •  1  >  )  ) 

T  ...  ;  1  •  J .  •  J u  ;  i 

v>  .  u  l  (  £  C  .  -1  * .  •  j  i ^  ) 

l  „  -  •  ~  ■  V-  r  )  Jw  i  (  f  *  i3  )  i  T  I  M 

t  .  L,  e  |)  !  ^  f  s  !I 

jl  .  i  t j  g: 

t  .  )  ♦  T  I  .  /  ■ )  uC  ) 


WRITE  (6,97) 

49  CONTINUE 

I F  ( MOD (  N  W  R  I  I  , 6 ) .  E  u  •  0  )  G^  To  80 
K  A  =  6  —  M  0  D  {  N  W  R  I  T  »  6  ) 

DO  58  1=1, KA 
DO  57  11  =  1 >  NCARDS 
57  WRITE  (2,95) 

53  CONTINUE 
8  0  STOP 
END 


$  I E  7TC  SO  5 1 

SUBROUTINE  FINDER  ( I , NCARDS »NP TS , I NkUT ) 

COMMON  A (8,2000) 

INTEGER  A 

4  FORMAT  (IX, 1814, 7X) 

DO  2  0  IA=1, NCARDS 
KB  =  18*1  A  -  9*( IA-1  ) 

KA  =  KB  -  17 

IF  (  I A  *  NE  * NCARDS )  GO  TO  11 
KB  =  KB  -  ( 18*NCARDS-NPTS ) 

11  READ  (INPUT, 4)  ( A (  I  >  IB )  » I B  =  KA *  KB ) 

IF  ( IA.EQ. NCARDS )  GO  TO  13 

KG  =  9 *  I A - 8 
KH  =  9* I A 
DO  12  I C  =  K  G , K  H 
KC  =  2* I C  -  KG 
KD  =  KC  +  1 

12  CALL  CuDER ( A ( I , IC )  , A  (  I  ,KC )  , A (  I  ,KD  )  ,2  ) 

GO  TO  20 

13  KE  =  9*NC ARDS-8 
KF  =  NPTS/2 

DO  14  I C  =  K  E  >  K F 
KC  =  2 *  I C  -  KE 
KD  =  KC  +  1 

14  CALL  CODlR  (  A  (  I  ,  I  C  )  ,•  A  (  I  ,  KC  )  ?  A  (  I  ,  KD  )  ,  2  ) 

2  0  CONTINUE 

RETURN 

END 

SI8FTC  SUB2 

SUBROUTINE  CODER ( I N , OUT  1 , OUT 2 , uP ) 

INTEGER  IN  ,0UT1 ,OUT2 ,OP 

C  OP  =  1  DECOMPRESS  •  . .  IN......OUT1  AND  OUT2 

C  OP  =  2  COMPRESS....  OUTl  AND  OUT2  . 

IF  (OP.EQ.2)  GO  TO  1 
OUT2  =  IN  -  IN/1G00Q*10000 
IF  (OUT2. GE. 90^0 )  OUT2  =  9000-OUT2 
OUTl  =  IN/ 10000 

IF  (OUTl. GE. 9000)  OUTl  =  9000-OUT1 
GO  TO  2 

1  I F ( OUT2 . LT . 0 )  OUT  2  =  -OUT2  +  9000 
IF  (OUT1.LT.O)  OUTl  =  -OUTl  +  9000 
IN  =  OUT2  +  OUT  1*10000 

2  RETURN 
END 

SENTRY  GE.NSTK 
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A3 *5  The  Synthetic  Seismogram  Program 

This  program  generates  a  synthetic  seismogram  given 
depth,  Z(J ),  and  velocity,  V(j),  values  for  a  layered  earth 
and  some  form  of  an  input  pulse.  For  the  input  pulse,  this 
program  either  calculates  the  Ricker  wavelet  given  by 

(4.2)  G(x)  =  -  |>5(x)  +  *6(x)] 


with  x  =  t/T  or  time  values  of  any  other  pulse  may  be  read 
in  as  data.  The  reflection  coefficient  Rq  is  computed  by 
an  iterative  procedure  from 


1-R 


v  . 

j 


(4.1) 


b.  np.  ,  +  (b.-b.  )  -b  .(3  .  (y-'-T.)  -p  .  ln- 

J-l  J-l  J  J-l  J  J  Vl+R-  j-1  v 

R  (od)  =  - J —  e  -1 

J*"1  b.  -.6.  ,  -  (b  .  -b  .  )+b.p.(j~~.^J) 


j-1" J-1 


j-l/  j  j  T+K 


J 


and  the  synthetic  seismogram  is  then  generated  from  the 
Fourier  integral 


(4.4) 


S(t) 
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Control  cards  required  for  use  in  this  program  are 
listed  below  (symbols  are  defined  on  comment  cards). 
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1.  NTSETS,  TINT,  MFREQ,  INTF,  INDEX  in 
(5X,  1515)  format. 

2.  (TE(l),  1=1, NTSETS)  in  (5X,  1515)  format. 

3.  (TST(l),  1=1,  NTSETS )  in  (5X,  1515)  format. 

4.  (AMP (i ) ,  1=1,  NTSETS)  in  (5(lX,Fl4.5) )  format. 

5.  DELAY  in  (lX,Fl4.5)  format. 


Cards  6  and  7  refer  to  SUB1,  subroutine  COEFF. 

6.  N  is  (1015)  format. 

7.  (V(J),  Z(J))  in  (IX,  El4 .7,  IX,  El4 .7)  format; 
one  card  is  required  for  each  set  of  values. 


The  next  cards  are  required  for  SUB2,  subroutine  PULSE. 

If  INDEX  =  1: 

8.  (FMT(J),  J=l, 12)  in  (lX,12A6)  format. 

9.  NS AMP  in  (5X,  1515)  format. 

10.  G(j),  the  pulse  amplitudes  in  the  time  domain 
in  format  FMT(j). 

11.  ITINT  in  (5X,  1515)  format. 

If  INDEX  =  2: 

8.  ITEND,  ITINT  in  (5X,  1515)  format. 

9.  PERIOD,  PULMID  in  (2(lX,  F9-5))  format. 
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£  JOB  7  3  5  vj  5  5 

£  I  J  J  0  3  R iv1  Civ 
£  It  rTC  SYNTH  NODcCK 

D  I  M E N S  ION  F W R  (  X 70  0  )  ?  F  W I  (  1  7 0 0  )  ,  T  I  ME  (  1 3 00  )  *  S  (  1 3 00  )  ?  Aivim  (  1  q  ) 
COMPLEX  RZ ( 1 70  0 )  »  FW ( 1700 ) 

COMMON  RZ , FW , OMEGA (1700) 

INTEGER  TINT?  TEND ? TSTART ? TE ( 10 )» TST ( 10 ) 

1  FORMAT  (  5  X  ?  1  5  I  5  ) 

3  FORMAT  ( 1 HL  >  4 ( 8  X , 4HT I  ME  ? 1 1 X  ?  6HS I GNA  L  ?  3 X )  ) 

3  FORMAT  ( 6 hi  END  ) 

6  FORMAT  ( 5 (  IX  9  F14 .3 )  ) 

7  FORMAT  ( 2 ( IX  ?  E14 • 7  )  ) 

3  FORMAT  (1X?4(2X?E14.7,2X?U4.7)) 

9  FORMAT  ( 1 H 1  9  34HAMP  L I  TUDE  MODIFIED  dY  A  FACTOR  OF  ?F14.5?13H  AT  Ti 
IE  T  =  9  F9 . 5 ) 

C  TINT  IS  (TIME  INTERVAL) *1000  FOR  CALCULATING  S  VALUES 

C  TST (  I  )  ARE  THE  TIMES  AT  WHICH  DIFFERENT  TIME  SECTIONS  BEGIN 

C  AMP (  I  )  ARE  AMPLITUDE  MODIFYING  FACTORS  INSTITUTED  AT  TIME  TST ( I ) 

C  TE (  I  )  ARE  THE  TIMES  AT  WHICH  SECTION  I  OF  PLOT  ENDS 

C  MFREQ  • •  MAX*  FREQ*  DESIRED  FoR  CALCULATIONS  up  FREQUENCY  SPECTRA 

C  INTF  «»  FREQUENCY  INTERVAL  FoR  WHICH  WE  WISH  TO  CALCULATE  SFcCTRA 

C  INDEX  =  1  READ  IN  PULSE  VALUES  (TIME  DOMAIN) 

C  INDEX  =  2  CALCULATE  GRAM—CHARL I ER  PULSE 

C  NTSETS  NUMBER  OF"  1  IMl  SECTIONS  Lou  divided  IimTO 

C  D  E  L  A  Y  *  *  *  *  D  E  L  A Y  TIME  To  GIVE  PROPER  START  FOR  AUIOPLOT  ^UTHJT 

READ  (  5  9 1  )  N  T  S  E  T  S  9  TINT?  MFREQ?  I  u  T  F  ?  i  N  D  E  X 
READ  (591)  ( TE ( I )  9  1  =  1  , NTSETS) 

READ  (591)  ( TST ( I )  9  I  =  1 » NTSETS ) 

READ  (5  9  6)  (AMP( I )  9  1  =  1  , NTSETS ) 

READ  (  5  ?  6  )  DELAY 

NTOT  =  MFREQ* INTF 

CALL  COEFF ( MFREQ  9  I NTF  9  NTOT ) 

CALL  PULSE (MFREQ, INTF 9 NTOT ? INDEX ) 

DO  11  J= 1 , NTOT 
F  iAi  (  J  )  =  R  Z  (  J  )  *  F  W  (  J  ) 

F  W  R  (  J  )  =  R  E  A  L  (  F  W  (  J  )  ) 

11  F  W I ( J )  =  A I MAG ( Fw ( j )  ) 

NF  =  NTOT-1 

DO  3 j  I A  = 1  ,  NTSETS 

TSTART  =  TST ( IA ) 

AMPMOD  =  AMP(IA) 

TEND  =  TE (  I A ) 

TMOD  =  FLOAT  (TSTART ) *  0 . 0  0 1 
WRITE  (6,9)  AMPMOD , TMOD 
WRITE  (6,3) 

DO  24  I=TSTART ?TEND, TINT 

T  =  FLOAT ( I ) *0. u01 

KB  =  (I  -  TSTART) /TINT  +  1 

IF  ( TSTART * EQ. 1 )  TIME (KB)  =  T+DELAY 

IF  ( T  S  T  AR  T • N  £  # 1 )  T I M E ( KB )  =  T 

A  =  0. 

B  —  0  • 

DO  25  J=1 , NF 

DFR  =  FWR(J-fl)  -  FWR(J) 

OF  I  =  FWI  ( J+l )  -  FWI  ( J ) 

WT 1  =  OMEGA ( J ) *T 
WT2  =  OMEGA ( J+l )*T 
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DELTC  =  C0SCWT2)  -  COS(WTl) 

DELTS  =  S  I  N  (  W  T  2  )  ~  SIN(WTl) 

A  =  A  +  DFR*UELTC 
2  5  B  -  8  +  DF I *DELT S 

DELTW  =  OMEGA ( 2 )  -  OMEGA ( 1 ) 

DEN  =  DELTW*T**2 
A  =  A/ DEN 
3  =  B/DEN 

S ( KB )  =  ( A~B)*AMPMOD*2. 

2  4  WRITE  (3,7)  T I  ME ( KB )  >  S ( KB ) 

IM  =  (TEND  -  TSTART ) /.TINT  +  1 

I F  ( MOD ( I M  , 4 ) . EQ . 0 )  I M4  =  I M / 4 

IF  ( MOD ( I M  » 4 ) • NE • 0 )  1 M 4  =  I M / 4  +  1 

DO  29  M=1 , 1M4 

M2  =  M  +  I M 4 

M3  =  M  +  2*1 M  4 

M4  =  M  +  3*1 M 4 

2  9  WRITE  (  6  9  8  )  T  I  ME  ( M  )  ,  S  ( M  )  ,  T  I  ME  (  M2  )  ,  S  (  M2  )  *  T  I  ME  ( I'M  3  )  ♦  S  (  M  3  )  ,  T  1  ME  (  i*i4  )  , 

IS ( M4 ) 

30  CONTINUE 

WRITE  (3,5) 

CALL  EXIT 
END 

6IBFTC  SJB1  NODECK 

SUBROUTINE  COEFF  (  MFREU  ,  I  jMTF  ,  NTOT  ) 

COMPLEX  RZ ( 1700 ) ,FW( 1700 ) 

COMMON  RZ , FW , OMEGA (1700) 

REAL  V ( 6  J 0 ) ,  Z ( 6 0 0  )  ,  S(600),  W,  VEL(6C0) 

COMPLEX  B ( 6C 0 )  »  R,  D,  F,  RZERC 
READ  (5,100)  N 
100  FORMAT  (1015) 

C  N  •••••••  THE  NUMBER  OF  SETS  OF  VELOCITY  AND  DEPTH  VALUES 

WRITE  (6,400) 

4  Oo  FORMAT  (  1H1 , 4X  ,  9HFRLUULi'>iCY  ,  8X  ,  22rlREFLECT  i  0  CotiFFIClcNI  ,8X,10riPHAS 
IE  OF  R,4X, 14HAMPLITUDE  OF  R) 

KK  =  N+l 

READ  (5,200)  ( V ( J )  , Z ( J )  » J  =  2 » KK ) 

2.o  FORMAT  ( 1 X  ,  E14 . 7 , IX , El 4 . 7 ) 

c  v(j),  z ( j )  the  velocity  and  depth  values 

V  (  1  )  =  V  (  2  ) 

Z  (  1  )  =  0  • 

5(1)  =  0.001 
S ( KK )  =  0.001 
DO  2  I  =  2  »  N 

S (  I )  =  ( V (  1  +  1  )-V (  I  )  )  /  (Z(  1+1  )-Z(  I  )  ) 

IF  ( ABS (S ( I )  ) .LT.0.001 )  S ( I  )  =  0.001 
2  CONTINUE 

DO  10  M =1 , N 
J=N-M+ 1 

IF  ( V(J+1 ) .EQ.V(J) )  GO  To  10 
VEL(M)  =  ALOG ( V ( J+l ) / V ( J ) ) 

10  CONTINUE 

DO  3  1=1, NTOT 

FREQ  =  FLOAT ( 1-1 ) /FLOAT ( INTF ) 

W  =  6. 2831854* FREQ 
IF  (FREQ.EQi^. )  GO  TO  7 
DO  4  J= 1 , KK 
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+  B(J)  =  CSQRT  (CMPLX  (  1  .~4.*(  W/S  (  J  )  )  **2  >0.  )  ) 

R  —  (  j  •  *  0  •  ) 

DO  5  M  =  1 * N 
J= N-M+ 1 
D=S ( J ) *B ( J ) 

E  =  S ( J  +  l  )  -  S ( J ) 

F=  S ( J+l )*B( J+l )*( 1. -R) / ( 1 .+R ) 

I  F  (  V(  J  +  l  )  .  EQ.  V  (  J  )  )  GO  TO  6 
R  =  ( ( D+E-F ) / ( D-E+F ) ) *CEXP ( -B ( J ) *VEL ( M ) ) 

GO  TO  5 

6  R  =  ( D+E-F )/( D-c+F ) 

5  CONTINUE 

RZERO  =  -R 
GO  TO  8 

7  RZERO  =  (  V  ( 1  )-V(i<K)  )  /  (V(KK)+V(1  )  ) 

8  RZ (  I  )  =  RZERO 

AM PL  =  CABS (RZERO) 

RR  =  RE AL ( RZERO ) 

R I  =  A  I  MAG (RZERO) 

PHASE  =  A TAN 2 ( R I » RR ) 

WRITE  (3>2oo)  FREQ  > AMPL 
3  WRITE  (6*300)  FREQ  >  RZERO , PHASE ,AMPL 
3 JO  FORMAT  (1X*5E16.7) 

WRITE  (3 >201) 

201  FORMAT  (6H  END  ) 

RETURN 

END 

S  I  BFTC  SUB2 

SUBROUTINE  P  UL  SE ( MFREQ  > I N TF , N  TOT , I NDE X ) 

COMPLEX  RZ ( 1 700 ) >FW ( 1700 ) 

COMMON  RZ , FW  >  OMEGA (17  0  0) 

DIMENSION  G ( 5 0 C )  >  FMT(12) 

1  FORMAT  ( 5X  > 1 5  I  5 ) 

2  FORMAT  ( 2 (  1X*F9.5  )  ) 

3  FORMAT  (6H  END  ) 

4  FORMAT  (lX>18F4.v>7X) 

3  FORMAT  ( 1  H 1  >  4  X  >  9  H  R  E  A  L  S  (  w  )  >  8  X  >  9  H  I  iv  i  A  G  S  (  w  )  >8X>9HAMPLI  TUDE>10X>5HPHA 
1 SE  >  1 1 X  >  9H FREQUENCY ) 

3  FORMAT  (  6 ( 2X>E14.7 )  ) 

7  FORMAT  ( 2 (  IX  >  E14 .7 )  ) 

430  FORMAT  ( 1 H 1  >  7X  >  4HT  I M  E'>  7'X  >  1 2  H  INPUT  SIGNAL) 

451  FORMAT (2X,F7. 3) 

452  FORMAT  (1X>12A6) 

GO  TO  (  1  j  >  11  )  >  INDEX 
10  READ  ( 5  >  45  2 )  (FMT (J)  >  J  = 1 > 1 2 ) 

READ  ( 5  > 1 )  NSAMP 

READ  ( 5  >  F M  T )  ( G ( J )  >  J  =  1 > NS A M P ) 

READ  (5*1)  ITINT 

C  NSAMP  THE  NUMBER  OF  AMPl  i  T  UDE  VALUES  >  G(J)  >  FOR  THE  PJLSt_ 

C  ITINT  FHL  DIGITIZING  INTERVAL  i  H  MSEC*  FOR  IimPUT  PULSi^. 

WRITE  (6>450) 

DO  16  J=l* NSAMP 

TIM  =  FLOAT (J*I TINT) *0*uul 

WRITE  (6  >  7 )  T  I  M  >  G ( J ) 

16  WRITE  (3*7)  T  I M  >  G ( J ) 

WRITE  (3*3) 

GO  TO  12 
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]l  A=(  1./SQRK  3. 1415927*2.  )  ) 

WRITE  (6*450) 

READ  (5,1)  I  TEND  , I T I  NT 
READ  (5*2)  PERIOD, PULMID 

C  ITEND  •  .  LENGTH  OF  TIME  FOR  WHICH  INPUT  PULSE  IS  TO  d  E  CALCULATED 

C  PERIOD  .  .  PARAMETER  CHARACTERISTIC  OF  PERIOD  OF  PULSo 

C  PULMID  ..  PARAMETER  FOR  PULSE  CALCULATION  PULM I D/PEK I OD  =  5 

DO  49  1  =  1  *  I  TEND  , ITINT 
TTT  =  FLOAT ( I - 1 > /1000. 

X  =  ( TTT-PULMI D ) /PER  1 OD 

FX= ( EXP ( -X*X/2. ) ) * ( ( ( ( ( (“X+l • ) *X  +  15 . ) *X“10. ) *X~45 . )*X+15- ) *X+15 . ) 
1A 

G  (  I  )  =  FX 

WRITE  (6,7)  TTT ,  FX 
49  WRITE  (3,7)  TTT  ,  FX 
WRITE  (3,3) 

NS AMP  =  ITEND 

12  SC A  LEI  =  0.001 
NL1  =  NSAMP- ITINT 
WRITE  (6,5) 

DO  20  I  F= 1  * NTOT 

FREQ  =  FLOAT ( I F  — 1 ) /FLOAT (  INTF ) 

W  =  6.2831854* FREQ 
SR  =  0. 

SI  =  J. 

I F  ( FREQ. EG. v. )  GO  TO  13 
DO  19  1=1 ,NL1 , ITINT 
IM  =  I  +  ITINT 
DELTG  =  G ( I M )  -  G ( I ) 

WT 1  =  W*FLOAT ( I )*SCALE1 
WT2  =  W*FLOAT ( I+ITINT )*SCALL1 
DELTC  =  COS ( WT 2 )  -  COS(WTl) 

D E L  T S  =  S I N ( WT 2  )  -  S I N ( W T 1  ) 

SR  =  SR  +  dELTG* DELTC 
19  SI  =  SI  +  DELTG*DELTS 
GO  TO  15 

13  DO  14  1=1 ,NL1, ITINT 
I M  =  I  +  ITINT 
DELTG  =  G ( IM )  -  G ( I ) 

TADD  =  (FLOAT ( I+ITINT )  +  FLOAT ( I ) ) *SCALE 1 

14  SR  =  SR  +  DELTG* TADD 
SI  =  SI/6.2831854 

SR  =  SR/6. 2831854 
F  W (  IF)  =  CMP  LX ( SR ,SI  ) 

GO  TO  21 

15  DEL  TT  =  SCALE1*FL0AT ( ITINT ) 

SR  =  SR/ ( DELTT*W**2*6 . 2831854 ) 

SI  =  — S I / ( D  E  L  T  T  *  w  **  2  *6 • 2  8  3 1 8  5  4 ) 

F  W (  IF)  =  CMP  LX (SR, SI  ) 

21  OMEGA ( IF)  =  W 

AMP  =  CABS ( FW (IF)) 

PHASE  =  AT  AN 2 (SI  *  SR ) 

WRITE  (3,7)  FREQ, AMP 
2  0  WRITE  (6,6)  FW ( I F ), AMP , PHASE  ,  FREQ 
WRITE  (3,3) 

RETURN 

END 


SENTRY 


SYNTH 


t4M. 


(  (  .  *K\  r  .  )  •  -  \  •  )  *  / 

T,  i  T  .  t  J  •  1  (  t  6  )  t. 

JO'I  T  j  H  . 

Tt  i  l  i  f  V  iT:  t  [»I  ^ 

. 

,  i  \(U1  U  «-!  TT  ,  =  X 

;  .  -/  (  .  -I  .-»«-))))  J  " ■  (  l  •  ‘A.  ■  :-)'*  >  ' 

x  *  tt  (V*d) 

(  «  ) 


T I  = 

.  tV  •  .  '  . 

(  .  .  • 

t  i  .  1 1 » i  j »  f :  =  i  i 
(,).=•  TJ 

T  » i  A  ,  J,  »  I-  I  I 

f  )  ■  :  i  :  )  T  ' 

.  \  o  s  i 

T  i  J  )  \  y 

(  j  f 

=  (  ) 

(  (  i  )  i  -  =  s 

54  .  r  _ 

<£*£)  TIP  >« 


